{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "### Advection-Dispersion Model for Solute Transport through a Column\n", "David J. Lampert, PhD, PE \n", "david.lampert@okstate.edu\n", "\n", "Laboratory columns are widely used to study the fate and transport of solutes in porous media. Tracer studies are often performed on columns to assess their flow properties, including their velocity and dispersion coefficient. Tracers are needed before columns can be used to study chemical processes in the column media to separate out the flow effects. This notebook provides a simulation tool for the flow a non-reactive tracer through a column. \n", "\n", "**Quick start instructions:** The notebook can be used by running the first code block to import functions and define parameters, then running the block to get the eigenvalues, then running the concentration calculator at the end.\n", "\n", "#### Model Problem\n", "Imagine that a 30-cm long column (a cylinder) has been set up with flow entering on one end and exiting on the other end. The velocity of the column is 10 cm/hr and the disperion coefficient is 100 cm$^{2}$/hr. Note that tracers are used to estimate these properties, but in this example, it is assumed that the parameter values are known. If the solute concentration in the influent for a tracer such as bromide is set to 100 mg/L, eventually the concentration leaving the column will have to reach 100 mg/L also. But what about beforehand? the concentration should start at zero and then increase eventually to the influent concentration level. To assess the flow properties, a mass transport model is needed. The diagram illustrates the model problem and the expected solution. The goal is to find a mathematical expression to estimate the breakthrough curve for this column.\n", "\n", "\n", "\n", "The parameters that affect the output from the column for this model include:\n", "\n", "|Parameter|Description|\n", "|:-|:----|\n", "|$C_{0}$|Solute influent concentration ($ML^{-3}$)| \n", "|$U$|Flow velocity in column ($LT^{-1}$)| \n", "|$D$|Dispersion coefficient ($L^{2} T^{-1}$)| \n", "|$L$|Length of column ($L$)| \n", "\n", "The code below defines the values of the function and imports key Python libraries that are used subsequently. " ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "# Define the parameter values in Python\n", "C0 = 100 # mg/L\n", "U = 10 # cm/hr\n", "D = 100 # cm2/hr\n", "L = 30 # cm\n", "n = 100 # number of terms to use in series solution\n", "\n", "# Import Numeric Python for array calculations and other math functions\n", "import numpy as np\n", "\n", "# Import Matplotlib for plotting the results\n", "from matplotlib import pyplot as plt, lines\n", "\n", "# Import the Scientific Python optimization library\n", "from scipy import optimize" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Governing Equations and Auxiliary Conditions\n", "The governing equation for advection and dispersion of a solute in one dimension is:\n", "\n", "$\\frac{dC}{dt}=D\\frac{d^{2}C}{dx^{2}}-U\\frac{dC}{dx}$ \n", "\n", "Two boundary conditions and one initial conditions are required to solve the equation. Assuming there is no tracer in the column at the start of the simulation, the iniital condition is zero concentration.\n", "\n", "$C(x,t=0)=0$\n", "\n", "The boundary condition for the influent is flux-matching, with advective flux into the colum matching the advective and dispersive fluxes at the start of the column.\n", "\n", "$UC_{0}=UC(x=0,t)-D\\frac{dc(x=0,t)}{dx}$\n", "\n", "The boundary condition at the effluent is normally assumed that the disperion flux is negligible, meaning the derivative is zero.\n", "\n", "$\\frac{dC(x=L,t)}{dx}=0$\n", "\n", "### Non-dimensionalization\n", "The equations can be non-dimensionalized in terms of dimensionless time $t^{*}$, distance $x^{*}$, and concentration $C^{*}$, which simplifies the mathematics. \n", "\n", "$t^{*}=\\frac{Dt}{L^{2}}$ $\\;\\;\\;\\;$ $x^{*}=\\frac{x}{L}$ $\\;\\;\\;\\;$ $C^{*}=\\frac{C}{C_{0}}$ \n", "\n", "In the dimensionless domain, the three parameters are reduced to just one parameter, defined as the Peclet number. \n", "\n", "$Pe$ = Ratio of advection to dispersion = $\\frac{UL}{D}$\n", "\n", "The governing equation, initial condition, and boundary conditions become:\n", "\n", "$\\frac{dC^{*}}{dt^{*}}=\\frac{d^{2}C^{*}}{dx^{*2}}-Pe\\frac{dC^{*}}{dx^{*}}$ \n", "\n", "$C^{*}(x^{*},t^{*}=0)=0$ \n", "\n", "$C^{*}(x^{*}=0,t^{*})-\\frac{1}{Pe}\\frac{dC^{*}(x^{*}=0,t^{*})}{dx^{*}}=1$ \n", "\n", "$\\frac{dC^{*}(x^{*}=1,t^{*})}{dx^{*}}=0$ \n", "\n", "### Analytical Solution\n", "The dimensionless governing equations and auxiliary conditions can be solved using [separation of variables](https://en.wikipedia.org/wiki/Separation_of_variables). The U.S. Geological Survey summarizes solutions for some important solute transport problems, including this one\n", "[here](https://pubs.usgs.gov/twri/twri3-b7/pdf/TWRI_3-B7.pdf). The solution is:\n", "\n", "$C^{*}(x^{*},t^{*})=1-2 \\textrm{ } Pe \\textrm{ } e^{ \\left(\\frac{Pe}{2}x^{*}-\\frac{Pe^{2}}{4} t^{*} \\right)} \\displaystyle\\sum_{n=1}^{\\infty} \n", "\\frac{\\beta_{i} \\left[ \\beta_{i} cos\\left(\\beta_{i} x^{*}) \\right) + \\frac{Pe}{2}sin \\left(\\beta_{i} x^{*}\\right) \\right]}\n", "{\\left[\\beta_{i}^{2} + \\frac{Pe^{2}}{4} + Pe \\right]\n", "\\left[\\beta_{i}^{2} + \\frac{Pe^{2}}{4} \\right]} \n", "e^{-\\beta_{i}^{2}t^{*}}$\n", "\n", "The values of $\\beta_{i}$ are the eigenvalues, and the corresponding terms in the infinite series are called eigenfunctions (see more description [here](http://tutorial.math.lamar.edu/Classes/DE/BVPEvals.aspx)). The eigenvalues are determined by finding the roots of the characteristic equation:\n", "\n", "$\\beta$ $cot \\beta$ $- \\frac{\\beta^{2}}{Pe} + \\frac{Pe}{4}=0$\n", "\n", "### Model application\n", "To apply the model, the eigenvalues must be estimated. The equation above has no exact solution, unlike some other characteristic equations that are commonly encountered in diffusion boundary value problems (e.g., $sin(\\beta L) = 0$ has $\\beta = \\frac{n \\pi}{L}$ with $n=1,2,3...$). \n", "\n", "The eigenvalues for a given column system with parameters $U$, $D$, and $C_{0}$ depend only on the Peclet number. The values for a given Peclet number can be determined by finding the roots of the function $F$: \n", "\n", "$F(Pe,\\beta)=\\beta \\cot \\beta - \\frac{\\beta^{2}}{Pe} + \\frac{Pe}{4}$ \n", "\n", "The function has a singularity at all integral multiples of $\\pi$, since: \n", "\n", "$\\cot \\beta = \\frac{\\cos \\beta}{\\sin \\beta}$ and $\\sin \\beta = 0$ at $\\beta = 0,\\pi, 2\\pi,... = n\\pi$. \n", "\n", "In between each singularity, the function has exactly one zero. The Python code below plots the value of the function across the first ten singularities. It also shows the first few roots." ] }, { "cell_type": "code", "execution_count": 225, "metadata": { "scrolled": false }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEhCAYAAACQrrywAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzsnXd4FNX+h98TeuhcQGkSiqLSQgBBBUTsCigiKiKKBVCuqNer/tSrEMEuXnsBLHhtqCgK2EVARVSSSFFQlKIYUOnSBIHz++PMbHaSTbLl7H6HMO/z7LO7szPnfPbMzPnM6UprTUBAQEBAAECatICAgICAAP8QmEJAQEBAQIjAFAICAgICQgSmEBAQEBAQIjCFgICAgIAQgSkEBAQEBIQ4oExBKZWtlNJKqQO+H65SaoibFkqpDAvhTXLCWpWwuP0MpVRGWFoOSWG8Sil1j1IqXym1z4k/M1XxF6NplaNjkqQOm4Sd22xpLamgzJiCUqqSUupfSql5SqktSqmdSqkflVLPKKWOkNYXD0nObNYBXzmvXTFoKu4GWe6E9U0iosLMJdKrcJwppYSHil0UpOW6FEo6E/g/oCGwwol/ezIiUkr1LOG8rArb9RtHx/Jk6AhIPuWlBdhAKVUbmAl0cDZtA34EmgCXAouBpSnSUlFrvTsVccWLo/Ed4B1bYWqtxwJjbYXn8FWh779aDt8KWuu1QFeBqFuHfW6vtU7YEKK8flfgNb+17getdb9ENQQIo7Xe71/AS4B2XvcDFcJ+Ox44wfmcHbbfMcB8YAeQB3QNO6YjxmTWYp4Ctzv7Xlgo3vA4JwFbgKnOby9gjGkrsBv4GXgEqFEojBOBD4HNwF/OMcOBIWHhh79mhx17PvClo2878AlwbNjvPcOOGwrMcuK4tlD4Gc7+hwFTgd+d/70G+Bg4pVBY4a9VzrGTwr872xRwBZDrpPM25/MxJZxLNxxdyjnvDixw/s83QLcwTdkR/n/PCOfN3a8q8Baw0knHXc55GANUdPaZXcz/HwJkhH8Pi6cN8Caw3rkGVgLjgGph+7jhzgb+CazCXDMzgINL+P8R9YT9fomT1jud//QlMCDs93DNNzr/fwfwUDHxhaflkBJ0rXL2mRS2rTXwmXOuvgf6FbNfdeC/Tjrtxtx/TwG1Ilwfq4ABTnjbgU+BVs4+A5199gAHhR17g7N9C1CF2O/zmK4rZ9vBwNNAPgX5wD1ApbB9jgI+cq6TXcBqzANbJ5H8VCJSq38AagJ/OydjAaBK2Dc77MRtdy4o99hVQHlnv3OAvc62PGBj2HFnRLgI3AtqMTDZ+W2bc5IXYIrS7r6vhx0/ANjnbN/pHL/JufDPwGR27nHLMTf2E86x/w777SfgF+fzbuDoCBfvLszT3RLgaiKbQq7zfaPz+Vfn+61AlhO/e8yvznfXBCe56Rj2/x4N238j8C0mYygpU3HD0SXscxAm43TTbQnmRo/XFOo6339z0nx12D73O/s8EZYe2vnvXzrnKSNs+xBn/yPCNG5zNO51vs8F0pz9Zoedt53AsrCwXiohDSLqcX67NWz7Lxhzd78Pd/YJ17zLSb/FwLhi4gtPy5LO3yrCMnugMgXX5t9OOuxwroPw/SpScP3tAhY66aad7RUKXR9/O2m2lIJ7aG5YnJucbVeFaXPDnxDnfR7rdfWPsPTY5vynXc736c4+aZj7UmMexvIw16GmkDmlLE+ViNTqH4DOYSfj0VL2zQ7bd6Sz7eqwbYc72xrgfcKojHly1MALES6CdUBjZ1s55z2zUNx3hF3MlZ1tK5xtK4GGzrbyQNsIN+6QsLDSw26Yu8Iurg+cbR9FuHhnhcVbjsim4GZi3cPiaozzBBbpwg/bPokwU3C0uzfr22Fx1waalXCOJoXFUfiV6exzu/N9H5DlbLu8sDaiv3krAkcW0vGCs8/qSNdPoX2LnCfgeQoePpo6264I26+Ps222831v2P9709n2W7TXc9i2qphM1033NKAC5klaA38428I1L8V5Gse5fiPEFZ6WhV8Phe23Cm9mf2nYfv2dbaeEbXP3u4iC+6O1s60p5mlfA4MiXB9uGv43bFsVZ9t45/vnzvdDw/Y5Js77PNbrapTzfQPQwNl2bNh+x2KMw/3eJCysljjXTapfZaGhWYV91jEc94LzviRs20HO+z7gAaXUGqXUHswTXEvnt4YRwnpDa/0rgNZ6r7PtBKXUt06Dtwb+42wvD9RTStUDmjnbJmmt1zjH79FaLy5Fe2vMzQ9wsxP+XuBkZ1uk+u3xWuu/CmkszHTnfaZS6gel1FuYp6n8UvREojMF5+a/YXFv0lqvjDKMrwq93Drzts77j1rrPOfzq3FodNkLXKiUWqaU2uWk54XOb5HOdzR0dt7naq1/dj6/HPZ7p0L7L9ZaL3A+u9dk/TjibY2pGgF4VWu9T2v9NzDF2VYPk9mG87zWejOUeG2E4zZqu69VJezbxnnfjTE7tNYfYJ7kw+nivJcHvnXOwSrMAwwUvaa3aK3d6zX8HnbTbJLzfoxSqgmmqhVgmdb6C+dzrPd5rLj/qQ6wxvlPn4f93lVrvQGY52pz8ozXMNXeayxoiJmy0ND8A+ZpojzQTSmltGO1JeHeBM6xLm4m9iKmrt99itoKHImp8yxHUX4L/6KUGoSpOwZTX7kaU0XR3NkWKYxYCDfC7zFF/3Ai/f/fImwrzEXANMyTUGvME92ZmAv0zJhVJojWOtHG2/B0KAeglKoZYb+bgJudzz9j0qox0IjEe+hF+6CyOeyze02qSDsmIe5oro1wxmqtJ8WqpZT70v2vf2OqUArze6HvkdIrFI7Wep5S6gegFXAeBaYwKWzfWO9zl2ivK/c/bQO+i/C7+x9OAC7AlByOBM7GVC23Aa4pQUdS2O9LClrrLcBrztcOwF1KqZDZKaV6KKV6xRismxlN1Fq3Bk7HnNhYj9+KqSrpgmlMDte9DlNtBHCxUupgR285pZT7dLUj7JCqYZ+/DfvtE0xxuKuTiQ4BRsegNZzumDaCK7TW3YG7nO3h6bczgp5IzKfg5rlWKVUJzM1jYVyEW5JqGdYvf0CE/f4I+9zCeY/UO8Y9X8u01hmYTggLI+wXOh9KqWj+P8CxSin3yfyCsN9zSjk+Xr6j4Bydr5RKU0pVwJT4wFR1/hzxyOTgnqtKSqk+AEqpUzDViOF87byXB64Nu567YaoLX4wj7ued939jMtt9wP/Cfo/3Po/2unL/k8a0D7j/6XjMQ+MbSimFud4maa0vdX53dceab9lBos7K9gtzgeVRUDf3J7AI09CrMRcZRK6D7Rl2XE9n21wK6nm/wxR13Uao2WHHeuoQw7YPDfstH1Pc3hC2LcPZL7yheYejeQMF9awq7D9sxRTV3baQG8PCW4tpIP0Dbz1tkf8WpnFIBD2/Ojp+cMLb7fw+N+w4N513YTI+t01jEmFtCs628IbmDc7/20GUDc0UNOa6r9ucfepT0P6xwzlHkRqay2MyQI252WdhMszC+90Ztm2lkw6h8xWmrW/YfqscTc2JrqH5O0puaA6/rrILx11MWkXcj6INzflh3yM1NBd7Poq5T9xOD6FX2H6r8F6DlcPOwW4KTKtwQ3MlCjpW7HP2W4qpMgy/N93rI/w6GxKmLSNse6OwNNfA+4X+U1z3OdFfV3XxdgBZhGmzcP97hhOWm2d9hzFRV3OxHQ2S+drvSwpg6qkxbvtvCtz5MMxN+TyFntKjYAgF3TfTMV04F8Vw/DOYxq/1mKLobEyjU2Hdr2PaAT7GFJtbYS7KL53fNcZgfsLUEx+FUx+stb4PGOTsWwPzfzdj/u/TMWgN51nMRfkPTPXROkx33/PD9rmagqe/Tk68xXE1MAJzs6djMtClmBsjGroUerUA0Fr/gXmqW4gpvu8B+hQ+WGu9B1N18I2zXx2gf4R47sKk22ZMWk7G9O4pzAxgIsYwmjqa0iMJ11ovBY7GdPHdhUmn1cADwCla632l/Pe40VrfgWngzcO0IdTGPFCcp7UebyGK5hQ9N8Vp+Qtzrj7HZHQVgcGYexOcUo3WehfGeP6LMZZDHe3fYTppfBurSK11Pqarp8ukQrsMIY77PNrrSmu9HlMaeRrzwHYE5vqaD9yCqRLbi+l2uwLTjnEY5qHkKUwX5ZSjHEcLCNjvCRtpfLvWOltSS0ABSqlDgZ+0W/xVqgcwx/l5uNZ6gpi4gCKUhYbmgIAAf3M/kKmUWoxpi+rubF9KfG0FAUmkTFQfBQQE+JpZmLafXpgeNqsw1UTHaq13lHBcgABB9VFAQEBAQIigpBAQEBAQEMKXbQp169bVGRkZohrWrFlDw4YN4ZdfYNMmaN9eVsfChVC7NhxyiJyG776DypWhRYvSD0qWjm3bYN8+OPxwOQ1798KGDZApt3TBphUrqL1pExx5JFSpUvoBSWBtfj4NfvsNGjeGgw4q/YAk8de331K5XDk4Qm6G/E3Ll1N782Zo3drcI0Lk5uau11rXSygQiX6wpb06duyopcnJyTEfRo7UumZNeR0HHaT1sGGyGtq21fqss0Q0hHSccILWxxwjq+Hqq7WuVUtMg9Za/3TPPVqD1osXi2nImzPHaBg3TkyD1lpv6tZNa+E8Y/kdd5i0+OEHUR1Ajg7GKSSZChXg77+lVUC5cuYJWZK0NJBug1JKXgP4Q0NAQBIITKEYOnVy5iurWFHUFEI60tJgbzRzlSVRg1KixtSpUydxUwhpEOb/brrJfBBMix7HHScWdziff/556TslmVtvvVVagjUCUygNt6Qg/WRYrpyYKYTwS0nBDwing6/KKdLXRIBVfNnQHIlff/2Vfv36kZeXx74UPa2q8AwoTc4/Qzr+9z/zktSQlyeaMasPP3QFWQszLS2NrKwspk6dSuPGjaMQ4RNjEiawgrLJflNS6NevH2effTY7d+4UbwgPXmXrtXPnTvr160e/fqUvLzx69GjzQctmiecOiDQpbGq55eabS98pBRx26KHSEuh/9tnSEqyx35hCXl4e//73v6lYsaK0lIAyRsWKFbn++uvJy4s0jb+X7OxsX5QUzjvvPPNB0JxuueUWsbjDOeywkuZkTA39+0eaZ3H/ZL8xhX379gWGEJA0KlasGFW1ZMOGzoJcwiWFyy6/XDR+gJYtW5a+Uwr46OOPxc/HP/8pMqFpUthvTCEgwA+sXbvWFyWFTZs3l75Tkvntd2cxNOEM+a9du0Tjh7DzIZwWNghMIQZ+++03zj//fFq0aMGRRx7J6aefzrJly5g9eza9e/dOqZa77rqr9J0icPnll7NkyZJif580aRJr1qyJev9wZs+eTc2aNcnMzAy9Pv7447h0FseqVat4+eWCpY5zcnK4+uqrrcYRFWXg5k+UIAUKKEtpEZhClGit6devHz179mT58uUsWbKEu+66i99/L7x0bOzs2bOn9J0KEY8p7N27l6effpojjzyy2H0Km0Jp+xeme/fuLFiwIPQ68cQTY9ZZEoVNoVOnTjzyyCNW4yiJrKwsX5QUmjdrZj4ImlOm0NQvhalZo4a0BJoJT8tjk8AUomTWrFlUqFCBK664IrQtMzOT7t3N1PDbtm3jnHPO4fDDD2fQoEFo52YdM2YMnTt3pk2bNgwbNiy0vWfPntxyyy0cd9xxPPzww0yfPp0uXbrQoUMHTjzxxJDZbNu2jUsuuYS2bdvSrl073njjDW666SZ27txJZmYmgwYNAuDFF1/kqKOOIjMzk+HDh7PXGdNQrVo1Ro0aRZcuXZg3bx49e/YkJyeHvXv3MmTIENq0aUPbtm158MEHmTJlCjk5OQwaNIjMzEx27twZ2h/g/fffJysri/bt23PCCSfElH533nknrVq14sQTT2TgwIGMGzculA5u+OvXr8ed82rVqlV0796drKwssrKy+OKLLwC46aab+Oyzz8jMzOTBBx/0lNI2btzIWWedRbt27ejatSuLFplFtLKzs7n00kvp2bMnzZs3T8hEcnNz4z7WJuMeeEBaAnPnzpWWAECPHj2kJcRdcvcj+804BQ/XXgsLFtgNMzMTHnqo2J+//fZbOnbsWOzv33zzDd999x0NGzbk2GOPZe7cuXTr1o2rrrqKUaPMSpyDBw9mxowZ9OljVo7cvHkzc+aYBag2bdrEl19+iVKKp59+mvvuu48HHniAsWPHUrNmTRYvXhzar3///jz22GMscNJg6dKlvPrqq8ydO5cKFSowYsQIXnrpJS666CK2b99OmzZtGDNmjEfvggULyM/P59tvvw1pqVWrFo899hjjxo0rGMXssG7dOoYOHcqnn35Ks2bN2LhxY8R0cDNslzfeeIPNmzczefJkvvnmG/bs2UNWVlaJaQlQv359PvroIypXrsyPP/7IwIEDycnJ4Z577mHcuHHMmDEDMFVWLqNHj6ZDhw689dZbfPLJJ1x00UWhNPr++++ZNWsWW7dupVWrVlx55ZVUqFChRA2RGDZsGBNq1BCvPnriiScYIaoArrrqKh4D8bRYuGgR7evWFdUwceJEhooqsEdQUrDEUUcdRePGjUlLSyMzM5NVq1YBpoTRpUsX2rZtyyeffMJ3330XOibUrRAzOO+UU06hbdu23H///aH9Pv74Y0/Phtq1axeJe+bMmeTm5tK5c2cyMzOZOXMmK1asAKBcuXIRu8s1b96cFStWMHLkSN5//31qlFIE//LLL+nRowfNnGqLOnXqRNyvcPVRixYt+Oyzz+jXrx/p6enUqFGDvn37lhgXwN9//83QoUNp27YtAwYMiKpd4/PPP2fw4MEA9OrViw0bNrBlyxYAzjjjDCpVqkTdunWpX79+3NV+EydO9EX10UeW22ri4dlJk6QlAPDLL79IS+CTWbOkJVhj/ywplPBEnyxat27NlClTiv29UqVKoc/lypVjz549/PXXX4wYMYKcnByaNGlCdnY2f/31V2i/qlWrhj6PHDmS6667jr59+zJ79mzTHx7TlqFKyYS01lx88cXcfffdRX6rXLky5cqVK7K9du3aLFy4kA8++IDHH3+c1157jWeffbbEOErTURLFHVu+fPlQV9DwtHnwwQc56KCDWLhwIfv27aNyFNMR6whPrG68kc5PQvilodkvOgLKDEFJIUp69erFrl27zJOiw/z580PVP5FwM7m6deuybdu2Ek1ly5YtNGrUCIDnn38+tP3kk0/mscceC33ftGkTABUqVOBvZ6K+E044gSlTpvDHH38Apm79559/LvH/rF+/nn379tG/f3/Gjh0bGrhVvXp1tm7dWmT/o48+mjlz5rBy5cpQHNHSo0cPpk6dys6dO9m6dSvTp08P/ZaRkRGqpw9Pny1bttCgQQPS0tJ44YUXQm0kxelz43nppZcAU61Ut27dUktAceGDkkJgBQHJIjCFKFFKMXXqVD766CNatGhB69atyc7OLhjMFIFatWqFqkDOOussOnfuXOy+2dnZDBgwgO7du1M3rH701ltvZdOmTbRp04b27dszyymmDhs2jHbt2jFo0CCOPPJI7rjjDk4++WTatWvHSSedZPrTl0B+fj49e/YkMzOTIUOGhEoZQ4YM4Yorrgg1NLvUq1ePCRMmcPbZZ9O+fXtP1Vc4bpuC+5oyZQpZWVmcd955ZGZm0r9//1DjPMD111/Pk08+yTHHHMP69etD20eMGMHzzz9P165dWbZsWahU1a5dO8qXL0/79u158MEHi6RhTk4O7dq146abbvKYqy3y8/PNB+En9GefeUY0foDlP/0kLQGAk048Ufx8PPH446LxW0V63plIr0iL7BipAWWB0aNH6/vvv19aRhGiucamTZum9Y03al25cgoUFc/Xt9yiNWj9zTdiGt55/XWj4Z57xDRorfXao47SukMHUQ05//63SYulS0V1ECyyExCQWkKN5MJPpne6XSAFdfT3waR8AF9//bW0BO73QRdhW1hvaFZK7dVaF23ZDAhwcBvR91uCNoWAMkwySgryd0xAQLIJev0ElFFiNgWlVH2l1BilVHHrzwV3S0CZZfz48b4oKYy48kppCTz26KPmg7BBtm/XTjR+gKE+mLXWFvGUFF4HfgbOB1BKtVFKjbOqKiDApwwbNsx8EM4ITznlFHEdfpi+G6Bp06bSEmKe9sXPxGMKVbTWzwB/A2itvwV6WVUVEOBTlFK+KCn0PessaQlUrlJFWgIA08LGvUhx/sCB0hKsEY8p/K6Uaoy3mqjUq0MpNbG0ffYH7rzzTlq3bk27du3IzMzkq6++iml66VgInywuWkaNGhWarvqhhx5ix44dod9OP/10NvtgHv6AgAD/Ek/vo38Bk4D6SqmBwMnA0rDflVLqtULHKOBUpVRNAK31uXHEK868efOYMWMGeXl5VKpUifXr17N7926efvppaWmAmRo7fOK7hx56iAsvvJD09HQA3n33XSlpZY+godlfBOfDGjGXFLTWPwGnA9cBRwI5wIWFdssA9gBPAI87r21hn/dL1q5dS926dUPz6NStW5eGDRt6nuirVavGf/7zH9q3b0/Xrl1DE68tX76crl270rlzZ0aNGkW1atUAiizQc9VVVzEpwkRjV155JZ06daJ169YFi8djpokYM2YM3bp14/XXX2fIkCFMmTKFRx55hDVr1nD88cdz/PHHh/Z1Rw1Hmmo70nTaAV569+7ti+qjzu4stoKZ4emnnSauAaD+QQeJxg+Q1aGDtARrRG0KSqlTlFI5SqlvgNuAt7TWt2mtH9da7wjbVWutjwLmAv8B/tJazwZ2aq3naK2LnyzI55x88smsXr2aww47jBEjRkSc92j79u107dqVhQsX0qNHj9BcSddccw3XXHMN8+fPL3FqjOK48847ycnJYdGiRcyZMye0VgCYSe8+//xzzj///NC2q6++moYNGzJr1qzQ1Bgu4VNtL1iwgHLlyvHSSy95ptNevHgxl1xyScw6yzqheZuEM8LbnOnYJXlz6lRpCQB07dJFWgI33nij+VAGSiyxlBQewZjBYKAekF3Szlrrx4ELgCuVUs8BsU9eXwLZ2dkopUKv3NxccnNzPdvcQVINGzYMbXPn8R82bJhn3/DVxoqjWrVq5ObmMmHCBOrVq8d5551X5Km+YsWKoSf/jh07hqbQnjdvHgOcEaAXXHBBzP/3tddeIysriw4dOvDdd9952jCKm4eoOIqbajvW6bQPRPr06eOLksKYsWOlJdCvXz9pCQB8+dVX0hK49/77pSVYI5Y2hR1a6/cAlFIjgC9KO0BrvQ4YopTqDnwfn8TIZGdnRxwZqyM4daQMf8KECUyYMCHmeMuVK0fPnj3p2bMnbdu2LTLpWoUKFULTNUczRXP41NHgnT7aZeXKlYwbN4758+dTu3ZthgwZUuwU3NGgS5hqO5bptA9EZsyYAe3biz8Rzp8/XzR+gHffe09aAoCpom3QQFSDO8twWSCWkkI9pdQFSqkOQGWgYrQHaq0/01rfG7M6n/HDDz/w448/hr4vWLAg6j7SXbt25Y033gBg8uTJoe1NmzZlyZIl7Nq1iy1btjBz5swix/75559UrVqVmjVr8vvvv/NelDdjcdNMFzfVdnHTaQcUwgclhRB+qK7wg4YAa8RSUrgP6AmMBI4AKiul3gQWAou01v6oYEwi27ZtY+TIkWzevJny5cvTsmVLJkyYwDnnnFPqsW5PoAceeIAzzjiDmjVrAtCkSRPOPfdc2rVrx6GHHkqHCA1W7du3p0OHDrRu3ZrmzZtz7LHHRqV32LBhnHbaaTRo0MDTrhA+1fa+ffuoUKECjz/+OFWqVOGSSy4JlVwilSQCHIQzQj9kw37QEJAE4p1eFdPDqDdwM/Bi2PZ9iU7dWhanzt6+fbvet2+f1lrrV155Rfft21dYUUBhor7GbrtNa6WSK6Y0pk83UzXPny+nYdcuo+HOO+U0aK11375aZ2bKapg82aTFkiWiMrAwdXbcs6RqrVcBq4AZhbYH03FHIDc3l6uuugqtNbVq1Qrq6vdTJkyYwDAQLym8//77nCqqAJ5++mn8MNHFqp9/JkNYw8cff8yJwhpsEWTgKaJ79+4sXLiQRYsW8emnn9KyZUtpSQFxMHz4cF+0KTzmrvQlaE4j/vlPsbjDWbBwobQEJvpkAKsN9htTSEtLY/fu3dIyAsoou3fvJi1tv7kd/FWfHzQ0lyms3gVKqWOVUpVshumSlZXFuHHjAmMIsM7u3bsZN24cWVlZ0R3gg5KCHwisIAJlwCBtPxq9BzSyHCYAU6dOZerUqVSpUsUz6Cx4Ba9EX1WqVAldX6Uxbdq0ZFzeMTPqttukJTD1zTelJQD+GNF8gzuiuQxg2xSS9gjVuHFj5s+fz969exNqWY/2lZ+fX/B9yBB0kyYpibdYHQ88gAb0li1yGoYNQx90kEg6hHRceim6cWOr4e7du5f58+fTuHHjUq9Dd0Q8IPpUGGqTEtQQdckqydSsVUtaAs2bNZOWYI39pxI1xTRqFFbgqVwZdu6U1VHOWfa6lBHSSdWQliaaCYV0SGtQ8tVHgy++WFoCTTMyzAfB8wHwwQcfiMYPcIUPVsKzRdxdUgGUUoVn5aoIXK2U2uhu0FqPYX+nShUxUwhR3jlVAqYQIi0NwqbkEEEp8UwohNa+MAgpfHIWAiyTkCkAhctMCmgM1HS+l43rxjUFyUzANYW9e2Xih8AUwjX4Bem0CChzJFR9pLW+JPwF7AJuDNt2qR2ZqWfo0KEFX6pUMZnh33/L6RAsKYQ0KCVqCkOHDhU3Bc91Iajj1FOlh67B5ZddJi0BgIymTcXN8YReZWdF4qBNoRg8M6i6a9EKVCGFdAiaQkiDcJvChAkTxJ/S/aABzHoZ0jz11FPSEgDI9MECN8OHD5eWYI3AFIrB08tE0BRCOgRNIaRBuPoopEPQmPzS+2jkyJFicbscddRR5oPwU3rhRaQkuOmmm6QlWMO2KdwFbCx1r/0Az7TRzhrHEqYQ0iFoCiENwqaQl5cnXn0U0iDMT8uXmw+CaZH7zTdicYezecsWaQmsWLlSWoI1Em1o9qC1LptzLQuWFEL4ofeRcJtCSINfGlcFdfgkBQIK45drMwGC6qNiaBC+kpOgKYR0CJpCSINwm0KDBg3ETSGkQZg6tWtLS6B0hHSQAAAgAElEQVTBwQdLSwCgcqWkzKwTE7V8MIDOFoEpFINnCU/XFHbskNMhaAohDcLVR2vWrBHPkD3XhaA5vfzyy2Jxu/z666/mg/DT8WmnnSYaPxDX0r5+JTCFYvCs/+yugSxgCiEdgiOaQxqETSGkQzATys7OFjcmgBdffNF8EEyL28f4Y1zq0u+tLv8eF6+99pq0BGtYNwWllODoKnvcfvvtBV/chmYBUwjpECwphDQItyncfvvt4tVHnutCkBdeeklagm9M4XsfmMLrU6ZIS7BGMkoK8o9RtnFLCtu3y2nwy4hm6Ya0oKE5IBLBubCG1d5HDqGzo5SqD1wF7NZa35GEuFKDYEkhhB96H7mmIDndhx9MwQfVRwEFBHZgl2S3KbwO/AycD6CUaqOUGpfkOK2Qk5NT8EWwpBDSIWgKIQ3uymRCmXJOTo54huy5LgTN6bFHHxXXMP/rr8U1ABx//PGi8QPcc8890hKskWxTqKK1fgb4G0Br/S2w/00SEpQUDK4pSI9VCEoKaB9o8EM6BNinVFNQSvVVSj2plDpFKZWmlBrivKLppPy7Uqox3hJe5bjVppBOnToVfKlQwfT+ESgphHQImkJIg5sJCJlCp06dxKuPPNeFoI6rrrpKLG6XTp07S0sAfDbNhfQDiwWiKSncC0wFhgJzgHOApsAMpVTfUo79FzAJqK+UGqiUeg6Q7yoQK0qZKiQ/NDT7oaQgeeEHbQoBPmT/t4IComlo3q21/lAp9QWwDqirtd6ulHoAmAkUu2it1vonpdTpwFlAWyAHeM6C7tRTrVpgCn6oPvKDKbj4QYcfNASUKaIpKXyklHoZOBr4p9Z6O4DWehtQN8L+yqlqylFKfQPcBryltb5Na/241lqwYj56Ro8e7d1QrRps2yanwx28JrCmQ0iDsCmMHj1a/CndDxoALrzwQmkJjB7lLLwobExHHH64aPwA5w4YIC3BGqWagtb6euA14Dzgn0qpn5RSM5RSHwOblVLNIxz2CMYMBgP1gGx7klODZ0QziJUUQjoExymENAi3KfhmRLOLoI6LLrpILG6XbJ8M5DviiCOkJXDuuedKS7BGVL2PtNZvaa0v11p3BFoBN2PaCmYBE5RSheeN3aG1fs/pbTSC/bDHUcOGDb0bhEoKIR0VKph3gZJCSINwm0LDhg3Fq49CGoQ5//zzpSUUvUeEePfdd8VLK54V+fZzYh68prXeCyx2Xi8Ws1tdpdQFwFLgB6Bi3AqFWLt2rXdD1aqwbp2cDtcUBNoUQhqEq4/Wrl0rbgqe60JQx4aNG8U1FLlHhNi5a5e0BDZt3iwtwRrJGNEMcD/QExgJHAFUVkq9CSwEFmmtpyYp3uRRrRqsWiUXv1t9JFBSCBE0NBdoEMZXzcvS5yPAKkkZvKa1fkRrPUxrfbTWuhZwOPAssBvon4w4bZOVleXdUK0abN0qp0OwpBDSINymkJWVJZ4hF7kuhDi0ZUtpCb5JCz+sZdC8WTNpCdawbgpa6yJhaq1Xaa1naK3v1lrLd5uIgtzcXO+G6tVF2hRCOgRLCiENwm0KIR2SS1CGXxeCOp588kmxuF2K3CNCnNBLvsnyvvvuk5ZgjWA9hWIYNmyYd0P16vDnnynPCEI6BBuaQxqEq4+GDRsmXn0U0iDMgw8+aD5Ip4UP8KynLsRTTz1lPpSBqrTAFIph4sSJ3g3Vq5vMMMVLcoZ0CA5eC2kQNoWJEyeKm4LnuhDUMePdd8XidilyjwixUrKtz+GjmTOlJVjDqikopXoopcTmNlq5ciWnnXYatWvXplGjRjz3nMXB09Wrm/co2xWsa1HKDGCLsaRgVUfQ0FygIaAA6fMRYBXbJYVZwCGWw4yac845h5NOOon169czceJE7rjD4hIOrin8+aeclgoVYi4pWNUh3NDs0eAHgszQP+cjOBfWsG0KYlfIokWL2LBhA9dddx3lnCkh6tWrF3d4+fn53g0xlBRsavHoKF8+ppKCLR0hDcINzSEdghlAfn6+LzLCVydPNh+k08IHnH766dISmDhhgrQEa5SZNoW5c+fSrVs39u3bR25uLtdddx1XXnll3OEV6VlRo4Z5j6KkYFOLR0eMJQVbOor0PhIqKeTm5opXH/ml99GyH38Ui9vFL72P/DBwbMWKFdISrFFmTGHBggV06tSJ448/nk6dOpGens7ZZ5/Nli1bOOqoo6hWrRrffvtt1OH17VtoVvCaNc17FKZQnJZ58+Zx9NFHc9xxxzFw4ED+juKp36MjxpJCcTp+//13jjnmGI477jh69epV6sjUkAZhU+jbt6+4KYQ0CHPrbbdJSyi4LoSrbr744gvR+AHuDlZe8x8LFiygc+fOzJo1i59++ok6depw4403kp6ezjvvvMM555yTWARuSWHLlri1NG3alE8++YQ5c+bQvHlz3n777dg0VKgQsylE0lG3bl0+//xz5syZw0UXXcQzzzwTXYB+aVPwS/2xX3RI4gODDLBLmTCFvXv3snTpUjp06EBaWhotWrTg2GOPBaBChQoJtS2EiLKkUJKWhg0bUqVKFQDKly9PWlqMyR9D9VFJOsqVKxeKe+vWrbRu3Tq6+P2yyI40ftDgEhhTgGXKhCn88MMP7Nixg/fee4+9e/eyYMECnnnmGS6++OK4wxw/frx3g2sKpZQUotGycuVK3nvvPXr37h2bjhiqj0rTsWDBArp06cJjjz1W6nQFIQ3C1Ufjx48vyJCFMkPP+RDMkP913XVicbsUuUeE8MN0G1cMH24+lAWT1lpbewH7gMMSDadjx446Fl588UXdtm1b3bRpU12tWjXdvn17/cYbb3j2ufjii/XixYtjCrcIFStqfeONCWnZsmWL7t69u/7+++9jj//ww7U+99yodo0mTbTW+tVXX9XDhw+PLv7Jk7UGrZcsiUW1XbKzjYa9e+U0PPyw0bBhg5yGTz4xGmbPltOgtdZpaVrfequshv79tW7dWlbDG2+Y87FwoagMIEcnmP/aLincCay3HGapLFiwgIEDB7Jq1Sq2bt3KggULOPvssxMKU0WqIqhVq9SSQkla9uzZw8CBA8nOzqZVq1ax64ihpFCSjl1hUw3XrFmT9PT06DQItykopcRLChGvCwGO98F8P0r4XLhMeeMN0fgBzu6/X8zzGRVWTUGbJTc32gwzGr755psSV186/fTT+fDDDxk6dCiTJk2KP6JataCU7m8laXnllVf46quvGDNmDD179uTVV1+NLf4YGppL0pGXl0ePHj04/vjjeeihh7jhhhuii99PbQp+KKYHGvzTviKdDmWIZK2nkFIWLlzI4SWs0/qurXliojCFkrQMHjyYwYMHxx9/DKZQko6jjz6aTz/9NPb4XVMQWBI0hB9MwQcZYZAFFhCkhV3KREPzunXrSjSFeIjYCByFKdjW4tERgynY1BHSIFxS6N27t7gpeM6HoDEd3bWrWNwu0XSUSAUNGjSQlkCnjh2lJVijTJhCMpg+fXrRjVGYQlJ1xDhOwboG4ZLC9OnTxU3Bo0GQu+++W1pCwXUhXHXTzelqLcktt9wiLcEacZuC8kuLW5Lo06dP0Y21a8OmTXI6KlYUMYWQBmf+JKmG5j59+oh3i/WcD8HM8OabbxbX0KdPH18Y5Odz50pL4K677pKWYI24TEEpNQDYqJTao5TKVUqdp5RKU0o9pJT6Sik1TimVYVVpipkxY0bRja4ppPBG9OioUAF2705Z3EU0CGfIM2bMEC8peDQI8sWXX0pLiHyPCFDaNC2pIMcn80DZIN6Swp3AY0BX4C3gOWAGcDEwB2gD5Cil2tgQaZNo5hsqljp1zJP69u32BMWCUPVRCL+spwDiVRa+0RAQYJl4TaER8IzWOkdrPRYYDpwC3KK1vlFrfSowHrC4oIEd2rdvz++//x7fwbVrm/c4q5D++OMPunTpEl/cYMUU9uzZw6GHHuoONoyNwBS8GvxAYEz+ogycj3hNYTkQ3rrzuvM+P2zbJExJwlccdNBBLFy4sNT9ImaadeqY943xDcVYv349f0a5SE9EHRUrJlx9tHbtWrZv3x7TIKyQBmFT0FqLm4LnfAhmAHPmzBGL2yWUFsIZ4YBEJ7u0wJtvviktwRrxmsI9wESl1N1Kqe6Y8Q6dgCVh+1QHqiSozzrt2rVj8eLFpe43IdKiGQmawo4dO0IT4kWLR4eFksLq1atp0qRJfBqETWHChAnipuDRIMi0adOkJfgmLVasWCFuTB9++KFo/DaJyxS01i8D52JKAp8Am4HJwDNKqRuUUr2BxwH5bgGFyMjI4Jdffil1v+HuBFfhuKawYUNccW/btq3UKSVK1GHBFH755ZeYTSGkQbhL6vDhw8VNIeJ1IcC4Bx6QluCbtMjJy5OWwJNPPSUtwRpxj2jWWs8AZiilqgDtgEygA9AfaIspJfymlHobWAQs0lq/Xlx4HpYvh+uug27d4LTTIMan65KIebrqcP7xD/MeZ0lhwYIFtGmTQNu7hS6pn376KZ07d47vYOEuqYC4KXgINASUQRIevKa13qm1/kprPV5rfYXWuium6qg1cD3wA6ZE8VjUgf71Fzz1FPTvD40awdixsHNnolLDNcd3oGsKcZYUvvjii9CaBnGRYJfUffv28fbbb9OvX7/4Aggamr0aBPGVFQTGVKZIyohmrfU+rfVSrfXLTm+kk7TWB0UdQOvWZjbSjz+GHj1g1Cho3x4sFBOVUlGZQsQ628qVIT09LlPQWjN37lyOOeaYmI7z6Eiw+ignJ4eaNWty2GGHxadB2BSmTZsmbgqe8yGYGd7jgxHNnvMhiC9GNLuDCcsA/p3mokIFOOEEeOstmDnTlBSOPRamTk0o2Gh73XQsbi6TunVhfeyzg69evZo9e/bQvHnzmI7z6EjQFN5++23OPPPMmI8LaRA2hY4dO4qbgkeDINFOvZ5Mir1HUkxtt6u4IC1atJCWYA3/mkI4vXqZUkJmJpxzDrz2WtxB/f3331EZQ6NGjSL/EKcpzJs3j2OOOSbm+fg9OipWNI28cWbKb731FmeddVbMx4U0CJtCo0aNxE3Bcz4ESwr93PVCBDUUe4+kmOk+GFl92eWXS0uwxv5hCgD16sFHH8Exx8CFF8Inn8QVTG5uLu3bt49fR5ymMGvWLI4++uj44wVTUoC4SgvLli1j06ZN8Tcyg7gpAOKm4NEgSFCL71PKQPtKwqaglDok0uR4ynBIouF7qFYNpk+HQw+FAQNg1aqYDtdaM2fOHHr06BG/hrp1Yd26mA7Jz8/ntddeY9CgQfHHC6akAHE1Nt9///1ccMEFifW+CtZT8OIHDX7AD+kgrMEHKWANGyWFlUC9CNvrOL/ZpVYtePttkzGdf35MT80rV65k3759tGzZstR9hw4dGvmHOEoKY8eO5bLLLouruO3R4ZpCjCWFWbNm8cEHHzBq1KiY4/doEO6SOnToUHFT8GgQpG+kWXxTjF/SolmM7XTJ4OSTTpKWYA0bpqCIbJTVgL8shF+Uli1hwgT46iu4556oD5szZw7HHXdcVPX6EUc0A9SvD3/+abrNRsFPP/3ElClTuOmmm6LWWayOOEoKO3bsYOjQoTzxxBPUqFEjMQ3BiGbv+RB8Or3xxhvFNRR7j6SYzp06SUtgxIgR0hKskch6Co8opR7BGMLd7nfn9TgwBVhgS2gRzj3XlBTGjoWlS6M6ZPbs2Rx33HFR7Vtsz4p6TqEoyiqk0aNHc8011/APd4xDjHh0xGEK2dnZdO7cOaFVsoLeR8VoEOTSyy6TluCb3kd+mGLi3//+t7QEayRSUmjrvBRwRNj3tkBLIA8YkqC+knn4YahaFUaOLDWT+O2335g+fXrkxXMikFfcmIj69c37H3+UGsaiRYuYOXMm1157bVRxlqojRlPIy8vj+eef5+GHH447fo8GYVPIy8sTNwXP+RB8Sv9h2TKxuF1CaSFcn78pxashRmL5ihXSEqyRyDQXxwMopZ4DrtFaxzb9pw3q1zclhZEjTQN0375FdtmzZw/ly5fn3nvvZfDgwYl3o3NNoZSSwsaNGzn//PMZO3Ys1atXTyxOlxhMYePGjQwZMoT777+f+q7mRAl6H3k1BARpUQaxMc3FJSKG4DJ8OLRqBTffXKRXzG+//UaNGjW49957ef7552Oq1y92MXA3gy1hTYYdO3bQu3dvzjjjjOIbrOPR4XZJLcUU8vPz6dGjByeffDKDBw9OKH6PBuHeRw0aNBA3BT8sEg9Q162OFDRHv6RFlcqVpSVQxwcD6GxhZZyCUuogpdQYpdQUpdTrSqnblVLRT2uRCBUqmNLCkiUwebLnp6+//hqAW2+9lYMPPjimkY9r1qyJ/MNBzt8qVH3kTnb3xx9/cN5559GiRQvuvffe6P9HNDqKKSls3ryZzz77DDDjEbp168ZFF13EuHHjYh4sV6IG4ZLCmjVrxE3Bcz4EM+S3fTB1drH3SIrpG6GGINU8++yz0hKsYWOcwrHAT8AFwE5Mj6NBwI9KqQRHa0VJ//7Qti3ceacnw5o3bx67du1iz549rFy5ko4dO7J69eqogszOzo78Q7VqZtbWQiWFq6++mqVLl5KVlcXu3bt59tlnExsTEElHMaZwySWXcMIJJ/DOO+9w3HHHceuttxb0TrFASINwl9Ts7GxxY8rOzvZFlckzzzwjLaH4eyTFfPvdd9ISmFzogXS/Rmud0AuYB0wA0sK2pTnbvognzI4dO+qYefllrUHrt94KberatavG9I7SgC5Xrpy+/vrrowrOJE0xNGum9YUXhr5+9tlnOj09XQM6LS1NP/jgg7Hrj0bH7NnmP86cGdr09ttv6/T0dJ2WlqYrVqyo33zzTWtxF9GwerWJf+JE63FErePZZ42GlSvlNDzzjNGwapWIBq21PtqUU7T+4AMxDYDWlStrfcMNYhq01noyaN2qlaiGM93zkZcnqgPI0Qnm6TaqjzKBB7TWoUc35/N/MesrpIYBA6BpU3AWH9Fah3pHpKen06VLFz788EPuu+++xOM6+GD47bdQPFdddRU7duwAzPTU//nPf1iyZElJIcRHocFrmzdvZsiQIezYsYN9+/ahlKKOuxBQMggamr0a/ID0aGIfpIUfRhP7QYMtbJjCFqBZhO3NMCuypYby5eHqq+Gzz2DBAn7++Wd2795Njx49mDNnDl9++SW9evWyUscebgr/+9//Qms+V61alRo1arBr1y5efvnlxOMpTKHqo4EDB7I5rDverl27uOyyy9ibrIbgwBS8CGrwwb8PKKPYMAV3Gc5BSqlmSqkMpdSFwETgFQvhR88ll5j6/iefJCMjg19//ZU5c+bQKY4Rjzk5OcX/ePDBsHYtYAbEtWzZkuuuu44nnniCmTNnsnnzZu644454/0XxOsJMYdmyZbz//vsopahevToNGjSgVatWZGRksNPigkQeDcKmkJOTI24KHg2CPOc2bAoaU4n3SArxwxQTD4wbJy3BGnGPUwjjRswAtmfDwvsbeBKIb26HeKldG847D15+GR54IHlT+x58sFloZ/dunnvuueTEEYkwUzjssMPYvXs35cuXt1P6iYZgQjwvQRWWwQ/nIsAeiTZKuC8gHTOauR2QnmBYOt7XMU6Dz5AEwijtdbkTR5MkxhHpleHEe3GK43VftZ34rxaKH9CDHA0tBTVc5GjIENTQxdFwiqAGQG8DfZ+whldALxXW0Nc5H5nCOrDQ0GzDDO4Eroiw/QpgbDxhxtX7yGXfPq0PPVTrnj3jD0OH9biJxLRpJum+/DKhOGLW8euvJt7x45Meb0QNmzeb+P/735TG79Hx4otGww8/yGmYNMloWL5cRIPWOmQK+t13xTQAWqenax1lj75k8QqI9z5yTSHofWQYDHwTYXsecJGF8GNDKbMIz+zZEOWYhJhp2NC8O+0KKaNSJfO+a1dq43UJGpq9GgQR/PcBZRwbplAfiDQR0HogNaOaCzNwoHlPYNnO0aNHF/+jawopGNHp0SFkCiENwqYwevRocVPwnA9BY7rcB7OklniPpJA2rVtLS+D8886TlmANG6bwC9A9wvYewK8Wwo+dQw+FDh3g9dfjDqLE0Zr165sMMgWm4NEhZAohDX4aTSyUIftBA4QtfCSoIXRdCDc0t2nTRlzDQPdBVFiHDWyYwnjgQaXUUKVUC+c1DHgAM6pZhgEDzCI8cVYhNXRLA5EoV870QEqBKXh0RDkhXtI0CJtCw4YNxTNkjwZBTj/jDGkJvkmLt95+W1oCQy65RFqCNWzMkvoAxhgeAZY5r4eBiVprC8OH4+Tss817nBOHrS2tvaBRI8jPjyvsuHUoZbqlprikENIg3CV17dq14qbgOR+CT4XrN2wQ11DqPZIidka5CmIy2bhpk7QEa1iZJVVrfTNQF+gKHA3U01qndoxCYVq1gsMOi9sUSiVFplCESpXkGpqFJ8QDxE3Bo0GQ/b+SIsCvWDEFAK31dq31fK3111rrbbbCTYjevU0vpG2xy8nKyip5h8aN4dfkN5kU0SFgCiENbmYoZApZWVnipuA5H4LGdHirVmJxu5R6j6SIWKbETxYtmjeXlmANa6bgS04/3dS/z54d86G5ubkl79C4MWzZEpfhJKSjYsWUtymENChlXkLVR7m5ueKm4AcNAC+88IK4htB1Idy4euopp4jGD/Dggw9KS7BG2TaFbt0gPR3efz/mQ4cNG1byDo0bm/ckVyEV0SFQUvBoSEsTKykMGzZMPEP2aBDkzrvukpbgm7T4ev58aQk8/vjj0hKsUbZNoVIl6NkTPvww5kMnTpxY8g6uKSS5CqmIDgFT8GgoV07MFCZOnChuCp60EHxCnvrWW+IaSr1HUsTy5culJfBBHHmMXynbpgBw4onw44/2Rze7ppCsUdPFIdnQDKakcKBPiOcDDUFDc0CyiNsUVMqm5kyQXr3M+6xZdsM9UE2hXLnAFPaTSz9l+GHAlh80lBHiMgWl1ABgo1Jqj1IqVyl1nlIqTSn1kFLqK6XUOKVUhlWl8dK2rZlS+9NPYzosv7S2gkqVzMjmJJtCER0CpuDRINimkJ+fXzBWQigT8KSFYEb0/nvviWvIz8/3hUGeedZZ0hIKptAvA+YUb0nhTuAxzLiEt4DngBnAxcAcoA2Qo5RqY0NkQqSlQffuMGdOTIeV2vsIoEmTpJtCER2VK0OKB+t4NAiWFDw9f4SMyS+9j5Z+/71Y3C5R3SMpYOPGjdISfNGuYYt4TaER8IzWOkdrPRYYDpwC3KK1vlFrfSpmlLOd5ccSpXt3+Okn+P33qA/p27dv6Ts1aQK//JKAsDh0CJQUPBoEG5r79u0rniF7NAhy7b/+ZT4IGlNU90gK+DTGWoBkMNbSSot+IF5TWA4cG/bdnXkuvG/YJExJQp5jjjHv8+bZDfeQQ4wppPLGFCgpeAgamgsIGpoDyiDxmsI9wESl1N1Kqe6YZTg7AUvC9qkOVElQnx2yssygr2SYwrZtZhBbqggams170NDsH/xg0AHWiMsUtNYvA+diSgKfAJuBycAzSqkblFK9gceBubaEJkTlypCZCV9+GfUh48ePL32nQw4x70msQiqiQ6Ck4NEg2NA8fvx4cVPwpIVgZnjrf/4jrsFzPgQ5qnNnaQn8c8QIaQnWiLtLqtZ6htb6eKAGpirpv8AWoD/wKtAFyFRKva2UGuv0WJKjSxfIzY36KbfUEc1QYAo//5yAsBh1SI9oFiwp+G5Es2CG3P+cc8TidonqHkkBLVu2lJbAqaeeKi3BGjamzt6ptf5Kaz1ea32F1rorpuqoNXA98AOmRPFYonElRKdOsH07RNlrI6phGE2bmvckmkIRHQIlBY8GwYZm5c69BGIZskeDIJkdOkhLiO4eSQEvv/KKeBVWH580utugfDIC1VrvA5Y6r5eTEUfMdOxo3vPywNbyffXrmyf3JJpCEaTbFIKG5gICDb7Q4IMUKFOU/WkuXA4/HKpUMaZgi7Q0U4WUSlOoXBn27JHLmIOGZl9o8E1G6JPSgm/wg0knyIFjCuXKQbt28M03Ue3eu3fv6MLNyEiqKRTR4a7TnMIqJI8GwYbm3r17i2fIHg2C9OjRw3wQzISivkeSTKOSls5NEZ190NhtiwPHFMD0QFq0KKobafr06dGF2bQprFqVmK5YdFRxevmmsArJo0GwpDB9+nRxU/CkhWCG/Oijj4rF7RL1PZJkevbsKS2BUaNGSUuwhhVTUEodq5SqVPiz72jfHjZtimq66z59+kQXZkYG/PEH7NyZmLZodVSubN5TWFLwaBBsaO7Tp4+4KfhBA8DIkSPF4naJ+h5JMrPjWETLNmPGjJGWYA1bJYX3MFNfFP7sL9q2Ne/fflvqrjNmzIguTLcHUpLGKhTRIWAKHg2CDc0zZswQz5A9GgSZ407tIGhMoetCuB49f80a0fjBHwv92MKWKahiPvsLt9fR4sX2wszIMO8rV9oLsyRcU0hSyaRUgobmAoKGZl8YZIBdDqw2hdq1oWFDWLrUXpiuKSSxXcGD26YgNf+RYEMz4A9T8FNG6AdzDChTHFimAHDEEfDdd6XupqO92Ro2hAoVklZSKKJDoPrIo0GwpKC1FjcFT1oIZsiLbZZ24yTqeyTJDLrgAnFznOGTRncbHJim8P33pV5EEyZMiC68tLSk9kAqokPAFDwaBBuaJ0yYIG4KftAA8Prrr5e+U5IJXRfCGfKPP/0kGj/A+++/Ly3BGgeeKbRqBVu3wtq1Je42fPjw6MPMyEiaKRTRIWAKHg2CDc3Dhw8Xz5A9GgS53e3tIpgh+yUtvvr6a2kJPPb449ISrHFgmgLAsmX2wmzWLHUNzdJtCkFDcwFBQ3NAYfxwXSbIgWcKhx1m3m2bwrp1Zm2FZOOH3kdBQ7O8Bhc/aAgoUyZtyxTuAjZG+Ow/mjQxU0X8+GOJu02bNi36MJs1M+9JmO6iiA4BU/BoEKw+mjZtmokfxHhySYoAABoVSURBVDLDadOm+aLK5LHHZCcdhhjvkSTS87jjpCUw6rbbpCVYw4opaK3v1lpvLvzZl6SlQfPmUMpC2x3dWVWjwTWFFSsSEBalDrf6KIWm4NEgWFLo2LFjQYYsqcFF8Cm9ta2ZfhMglBbCpZU6//iHaPzgjzUdbHHgVR8BtGgBpfRYaNQohkHZrikkoV2hiA4BU/BoECwpNGrUSLzqxg8aAHoef7y4Bk9aCPLmm29KS+Ciiy+WlmCNA9MUmjc3GbitG6pePUhPT01jsx/aFIKGZnEO7H8fkEwOXFPYtg3Wr7cTnlKp64GUlmbaRIKG5gNbg4sfNPiBIB2sYWuW1EpKqWZKqSOVUvVshJlUoqjuGTp0aGxhNm+elDaFiDqqVEmpKXg0CFYfDR06VDxD9mgQZMAA2SXPIey6EM6Q/VCff+opp0hLsEbcpqCUqq6UulIp9SmwBfgJ+Bb4TSm1Wik1USnlz5UnolhbOeoRzS5uScHyDRJRR+XKKTWFYERzIQ0ugpmhH6Zq9pwPQbp27SotwRdTmdsiLlNQSv0LWAVcCnwEnAlkAocBRwOjMes/f6SUel8pdagVtbaIwhRi6n0EBVVSGzYkICxKHSkuKXg0CJYUPL2PhDJkP2gA6N+/v7iGmO+RJPHOu+9KS+Caa64xH8pANVb5OI87BjhOa13cwgRfA88qpa4ALgOOA0oeGJBKatWCGjVKXAMhL9a1nMO7pdatm4C4KHSk2BQ8GgQbmvPy8sQzZI8GQb5dskRagkkLH3QH3bhxo7mnBfmxlC7u+xNxmYLWOqoKTa31LuCJeOJIOk2awOrV9sJr3ty8r1wJRx1lL9xIpNgUPAS9jwrwg4aAAMsk0qYg/7iUCE2alFhSaNCgQWzhJWmxnYg60tNTagoeDYKm0KBBA3FT8IMGgPr16olrCF0XwuZYxR27I0id2rWlJVgj3jaFAcAmpdQepVSuUuo8pVSaUuohpdRXSqlxSqkMq0pt06RJiWs1r4l1ib9q1cx4BcvFyIg60tNhxw6r8UStQdAU1qxZIz7NxZo1a3xRffTZ559LS/BNWpzjtq8I8uKLL0pLsEa8JYU7gUeBrsBbwHPADOBiYA7QBshRSrWxITIpNGoEf/wBu3dH/Dk7Ozv2MN1BcRaJqCPFJQWPBkFTyM7OFp/mwpMWgk/Ijz76qLiGuO6RJLBw0SLx0spLL70kGr9N4jWFRsAzWuscrfVYYDhwCnCL1vpGrfWpwHjgDks67dO4sXkvpkRw++23xx5mEsYqRNRRpUpKSwoeDYKmcPvtt4tX3Xg0CPKIDybEi+seSQILFy2SlsBLL78sLcEa8ZrCcuDYsO/uMlDzw7ZNwpQk/EnDhuY91mqikmjWzLRT7NljL8xIpLj6yIN0Q7Nw9ZEHP2jwA0E6lCniNYV7gIlKqbuVUt0xvZg6AeH95KoD8i1AxeE2kpWyAltMNG9uMkybvZoicSCbgnD1kUeDHzJDaQ0+KDUF2CUuU9BavwyciykJfAJsBiYDzyilblBK9QYeB+baEmqdUkwhJycn9jDdbqkWq5Ai6khxm4JHQ/nyyS8JlaRDOEP2aBBk6tSp0hLiu0eSwBmnny4tgYcfekhagjXiHbyG1noGMEMpVQVohxnR3AHoD7TFlBJ+U0q9DSwCFmmt5Vcbd6lXzzz12i4pQPInxqtSxTSQ79ljMulUIl1SCKqP/KchoIAycD4SnhBPa71Ta/2V1nq81voKrXVXTNVRa+B64AdMiUK+ZSyctDSoXx9+/z3iz506dYo9zMaNTSZtsaQQUUd6unlPURWSR4OgKXTq1Em8+sijQTADOLNfP7G4XeK6R5KAH6a5uPraa6UlWCMpj5la633AUufl32b5gw4q1hTiolw5M69SsksKrins3Gmm60gl0iUFH2TIfqg+8hVl4Ok4oIADcz0Fl4MPtmsKkJp1FapWNe/btyc3nki4s6RKZQRB9ZG/NAQGWeY4sE2hhOqj0aNHxxdms2ZWq48i6khx9ZFHQ7ly5l2g+mb06NHi1UceDYIZ8tVXXy0Wt0vc94hl2rdrJ26Ogy64QDR+myTVFJRShyil/Gs89erBunURL6i4R2s2a2bC3LYtMW0l6XBLCikyhSIjmkGkCskzolkoE/BoEMQPUzX7ZURz+8xMaQlceOGF0hKskewMexWwUCnVI8nxxEf9+qZePkI1TEN3cFusWO6BFFGHW1JIUfWRR4OgKTRs2FC8+siTFoIZchcfLCwT9z1imddel+/UOGjQIGkJ1ki2KVwKvAncn+R44sOdaXLduiI/rY23q2r4ugoWiKgjxSUFjwZBU1i7dq149ZFHgyB/RLhmU03ouhCuutkpNY18GBs3bZKWYI2kmoLWepLWerTWuksy44mbEkwhbqJY/zlhUtym4EHQFADx6iMPgQZfGGSAXfxb358K3FWjIiyhmZWVFV+YdeuaabQtmUJEHSnufeTRIGgKWVlZ4qbgBw0ArVu3FovbJe57xDL/qFNHWgItW7QwH6RN2gIJm4JSqr5S6g2l1Gal1B9KqcY2hKUEd9nM9euL/JSbmxtfmEpZ7ZYaUUeK2xQ8GgRNITc3V7xNITc31xdPxzNmzDAfBDOhuO8Ry/Tu3VtaAo/6YNZaW9goKTwB1AHOAaoBFQGUUo8ppf7PQvjJwzWFCCWFYcOGxR9uRgasWhX/8aXpSHFJwaNB0BSGDRsm3qbgSQvBDPmmm28Wi9sllBbCT8dfzJsnGj/Aww8/LC3BGjZM4QTgGq31x0B4TvE2cL6F8JNHzZrmyTNCSWHixInxh+uWFCzcLBF1VKliMkdL3V5j0iBoChMnThSvuvGDBoBXJk8Wi9vFkxaC/Pjjj9ISeO/996UlWMOGKewF/oqwfTnQ3EL4ySMtDWrXBts9B5o1Mxl2hBKIFdLSTBWS1IhmkGto9sOIZh9khCHKQB22FYJ0sIYNU3gHuCjC9hp4Sw7+pE4d2LjRbpip6IFUteqBaQp+WE/BRTAjCrLAAoK0sIsNU7gFuFwpdQegAK2USgdGAXkWwk8uxZhCfn5+/GFmZJh3C+0KxepIoSl4NAiaQn5+vnjVjR80AMz/+muxuF0Sukcscu6AAdISeOnFF6UlWMPG1Nn5mKmxOwPpGCPYjFmu88ZEw086xVQfJdSzwmJJoVgdKTQFjwZ3/QaBhXaC3kcFLF682HzwQ+8j4aqbDcmqpo0BP7Rr2CIhU1BKtVdKDQGOAPoATTFVSX2BVlrr/bak0Ldv3/jDrFHDhGvBFIrVUbVqyhqaPRoESwp9+/YVrz7ypIVgZnjJZZeJxe3iOR+CzPzkE2kJZN9+u7QEa8S9noJSahjwJKbKCOBHoJfWeroNYSmjVi3YvNl+uBkZ8PPP9sN1qVbtwG5TCBqaDUEDa4BlEikp3IgZo3AwpuroD+BeG6JSSs2axhRs31wWxypEpFq1lJUUPEibgh96H7kEDc0BhfHDdZkgiZhCU2Cc1voPrXUuMAQ424qqVFKrlsngCs0jNH78+MTCdU0hwYukWB0pNAWPBrdNQcAUxo8fL1595NEgmAHce889YnG7hK4L4Yzw6KOPFo0f/LG+hS0SMYVyQGh6Qq31cgClVINERaWUWrXMe6EqpIRGNINZlnPnzoQn2ytWRwpNIeKIZoGGZs+IZqGMyKNBkNBUzYIZsl/SotVhh0lL4PTTT5eWYI1Eex8NU0r1Ukq5M1LtBaokGGZqcU2hUA8klejF7nZLTbBdoVgdKTQFjwbBkoJSSrz6yJMWghlyoyZNxOJ2SfgescSk558XL62ccuqpovHbJBFTmA1cB3wMrFNKrQYqY4ziJKVUbQv6ko+78P2ff9oN1zWFZA1gq1bNVHmluhpFsKQAiFcfeTT4of7YDxqECVLALnH3PtJa9wJQSjUHOoa9Lsc0Qmul1Aqt9aE2hCaNmjXNu21TaNrUvCerB1K1auZ9+3aoXj05cURCsKQA+CND9sETcpARBiQLG4PXVmitX9da36S1PklrXRcz59F5gPw6eaVRTEkh4el4a9Y0VVMJ9kAqVodrCimoQvJoECwp9O7dW7z6yJMWgsZ00oknisXtEkoL4dJKk8bys/V3OeooaQnWiMsUlFLNSvpda71Kaz1Fa32LMshXgBZHMaYwfbqF4RZNmyZcUihWRwpNwaNBsKQwffp08eojjwZBnn/+efNBMEP2S1qc6AODHDNmjLQEa8RbUpinlHpGKVVsXzClVG2l1JXAEuDMOONJPsWYQp8+fRIP24IpFKvDrTLaujWh8GPWIFhS6NOnj3j1kSctBDPkiy6+WCxuFyv3iAU+/vhjaQmMGjVKWoI14jWFw4GNwDtKqXVKqfeVUs8ppZ5USk1WSi3CDGa7ELhWa+3fZYnczHXLFs/m0MpWiXDIIcYUEsg8itWRwpKCR4NgSSGkQymxDHnGjBnixgTwkZsRCmqwco9YYPWvv0pL4Ct3gsIy0PAflylorTdrrW8AGgFXAt8DtYBmwB7geaCD1vpYrfUHtsQmhbS05M0j1LSpeZIvZDhWSGFJwYN07yMwmbIfeh8J4quspwxkhIlSllIgkbmPGjkzpE5xXvsv1asnJ3N1eyD98kvBeAhbpLCk4EG69xGIlhQ8+EGDND4wSCA4FxZJpPfRL0qp+taUSBJhIJi2cZFZ6JZarI4UlhQ8GgRLCiEdaWlimYDW2hfVR7+tXSuuwco9YoFLLrlEWgIffuDvCpFYSMQUPI8ISqnF+90UFy4RSgoTJkxIPNxDDjHvv/wSdxDF6kihKXg0CJYUQjoEq48mTJjgi6fjF3ywqIuVe8QCP/zwg7QE3nn3XWkJ1rCx8ppLBvvbFBcuEUoKw4cPTzzc+vWhYsWETKFYHW71UQpMwaNBsKQQ0iFYfeRJC8En5etvuEFcg5V7xAJzv/hCWgIPP/ywtARr2DSF/Zdq1ZKTuaalQZMmyRnVXL48VKmS+obmoE3BF9VH/qi4cfBJNVKAHRI1hUuUUl2VUpUx1+n+eXUkq6EZTBXS6tXJCbtGDfvTc5SGH3ofCbYpAL6oPvINQVqUORKdEO//gC+AP4FqwD1KqZFKqWOVUlUt6EsNVasWWU9h2rRpdsJu0iQhUyhRRzLNrDgNgovshHQItil40kLQmF743//ENVi7RxLkxBNOkJbAmGA5zhInxBsN1AH2KaWWaa2PtCE0qVStWmRpy44dO9oJu0kTyM83T9blY0/uEnWkyBQ8Gtz/IFBSCOkQrD7q2LFj6qvsItA+M1Nagr17JEHq1q0Ly5aJajj0UH/P+xkLyZwQ73zgrYQVpoIIptCoUSM7YR9yiHmqdbsQxkiJOmrUSEkG5dEgWFII6RCsPvKkheBTett27cTidgmlhXCbwuRXXxWNH2DgBReYD2WgfSXukkJJaK1XAavYXwa1Va0Ku3fH/TRfIu5iKL/8UvDZFtWrm1JIKhEsKYTwy4hmP2QA0hr80qYgnA4+uBKsEfQ+AkhPN++FSgtWcI0gGfOzSDY0B72PRClLmVCiBGlhl8AUwJQUwGMKQ4cOtRN2gqZQoo4UmYJHg2BJIaRDsPrIkxaSs6QOHiyuwdo9kiC+WKP5tNOkJVgjMAWIaArWRmvWqGHGQcTZA6lEHSlqaPZocBe5CUY0i2bI/33wQbG4Xfwyorlbt27SEvjXv/4lLcEagSlARFOw1rNCKVNaiLOkUKKOGjXgr79Me0gS8WhQylQhHci9j3xQfdSrVy9pCQXnQ7g+/6233xaNH2DEiBHSEqwRmAIUtCns3BnalJeXZy/8xo3jNoUSdRSzQJBtimgoX16kpBDSIVh95EkLwcxw4aJF4hry8vJ8YZAbNmyQlsCPP/0kLcEagSmAmS4Cigxgs0YCplAirilITHUR9D4KprkIKJMEpgARSwoNGlic8LVRIzNOIY6MtEQdNWua92Qs4lOSBiFTCOkQrD5q0KCBL56ODz7oIPNB0Jis3iMJkF5Ffh7Of9SpIy3BGoEpQEFJIcwU1qxZYy/8xo3Nk+1vv8V8aIk6UlR9VESDkCmEdAhWH3nSQjBDXrJ0qVjcLqG0EG5TGDhwoGj8AJMnTzYfpMeNWCAwBSgoKYRVH2VnZ9sL3x35GcdAsxJ1pMgUimgQMoWQDsHqo+zsbF+UFO655x5pCb5Ji7y8PPHM+H8vvCAav00CU4CIJYXbbU5w1bixeY+jXaFEHW71UZJNoYgGIVMI6RCsPvKkhWBGdO9994lrsHqPJEDeggXSEgJTKHNEMAWrJFBSKBG3pJDkNoUiSDc0+2Xq7KChOaAMEpgCRKw+skrdumYFNtumkKKG5iJIm4Jfeh/5gTJQhx3gLwJTAKhQwTx9hpUUcnJy7IWvFDRsGJcplKijcmWjPcmmUESDkCmEdAhWH3nSQjBDnjV7tljcLqG0EDams848UzR+gCcef1xagjUCU+D/27vfGLmqOozj3yf9k25LTY2oYdsqEMgiIu5uFlLFbBqrpiCCGmPQQDTBmJhqQF8Q/7w3MRriOxOlGhOxiqUmWhKgRqqSSAuUtmxZCi1UWP61StGtmFDNzxfnzHWw292mnTnn0j6fZLN3ZrNzn71z5/zuObNzDqmRGRhInw7ul8HBk54++7ik1Fs403oKHj5qjzb1mqwnXBQ6Fi16XVEYGxvr7eMPDsJJ/JvrnDkKFIVjMlQqCk2OisNHY2NjrWgIV69enTYqFqaev0ZOUiumuVi3rnaEnnFR6Pi/otBzJzl8NKczsadQe+rsDr/RbKchF4WOEkVhehqOHOnt49YqCkePlt1nNw8f/Y8ztCcDtCfHKVC08I+QNA3srRzjbOCvlTNAO3K0IQO0I0cbMkA7crQhA7QjRxsyAAxFxNJTeYC+LMfZA3sjouqApaSHamdoS442ZGhLjjZkaEuONmRoS442ZOjkONXH8PCRmZk1XBTMzKzR1qLQhnX+2pAB2pGjDRmgHTnakAHakaMNGaAdOdqQAXqQo5VvNJuZWR1t7SmYmVkFLgpmZtZwUTAzs4aLgpmZNaoXBUlrJe2VtE/S1ytl+LGkg5Imauw/Z1gp6T5Jk5L2SLqpQoZFkrZL2pUzVF1aS9I8SY9I2lxp/wckPSppZy8+FHQKOZZJ2ijp8Xx+vK/w/ofyMeh8/UPSzSUzdGX5aj43JyRtkLSoQoab8v73lDwOx2unet6GRkS1L2AesB84H1gI7AIurpBjHBgFJioei3OA0by9FHii9LEABJyVtxcA24BVFY/J14CfA5sr7f8AcHatv78rx0+BL+TthcCyilnmAS8C76yw7+XA08BAvn0H8PnCGS4BJoDFpBkhfgdcWGjfx7RT/WhDa/cULgf2RcRTEfEa8AvgWkkXSDqUr9R2SnpZ0n5Jb+pHiIj4I/By930VMrwQETvy9jQwCSwvmSOSzox9C/JXlD4WAJJWAB8Fbuu6r3iOGXIVzZAfdxxYDxARr0XEKxWPxRpgf0T8pVKG+cCApPmkhvn5wjneBTwQEa9GxL+BPwCfKJFhpnaKPrShtYvCcuDZrttTwPKI2AfcD9wQEcPAbuDjEdHfFeq71Mwg6VxgBNhWOkcestkJHAS2RETxDNn3gVuAZuGECjkCuFfSw5K+WCnD+cAh4Cd5KO02SUsqnp/XARug/LGIiOeA7wHPAC8Af4+IewvnmADGJb1F0mLgKmBlxeej521o7aIw02olnU/TvZv0BABcRJ1ZU4tnkHQWcCdwc9cTWCxHRPwnn0QrgMslXVI6g6SrgYMR8fAMPy75nFwREaPAlcA6SeMVMswnDRn8ICJGgH8CnXHjouenpIXANcCvuu4ueV68GbgWOA8YBJZIur5kjoiYBL4DbAHuJg3XdBYXqdFm9bwNrV0UpoCVXbdXkLqDA8CiiDgsaSXwt9w1KqZGBkkLSAXh9ojYVCsHQES8AmwF1lbIcAVwjaQDpO7wByX9rHSOiHg+fz8I/JpUJEsfiylgKiK25dsbgdFK58WVwI6IeAmqnJsfAp6OiEMRcRTYBLy/wnmxPiJGI2KcNJzzZMU2q+dtaO2i8CBwoaTz8lXIdcBvgItJY+qQxvAmj/P7/VQ0gySRxo0nI+LWGjkkvVXSsrw9QHoRPl4yA0BEfCMiVkTEuaRz4vcRcX3JHJKWSFra2QY+QrrqKn0sXgSelTSU71oDPFY6R/YZ8tBRVjrDM8AqSYvz62VN3mfp1+rb8vd3AJ8kHZNabVbP29CqRSG/UfNl4B5S6DsiYg+v7/b8i3RldFG/ckjaAPwZGJI0JenG0hlIV8c3kK6KO//6d1XhHOcA90naTTrZtkTE5sIZZlMyx9uB+yXtArYDd0XE3YUzdHwFuD0/L8PAt0vnyOPnHyZdnXcUzZB7SxuBHcCjpPbrh6VzAHdKegz4LbAuIg6XyDBTO9WPNtQT4pmZWaP28JGZmbWIi4KZmTVcFMzMrOGiYGZmDRcFMzNruCiYmVnDRcHMzBouCmYnSNJnldZXeFXSk5I+XTuTWa+5KJidgDxJ33rgu6Q59X8J/EjSvKrBzHrMn2g2OwGStgLbI+KWfPsy0iJE88IvIjuNuKdgNoc8788HgLu67l4L7HJBsNPN/NoBzN4ALiVdQD2SZ4/9FPBN4Maqqcz6wEXBbG7DpHVwLwAeIi1scg9pjWCz04qHj8zmNkKarvkJYBVpquJVwK2z/ZLZG5F7CmZzGwY2RcQR0toK2/M62qsrZjLrC/cUzGaR/+X0PRy7ctWlwJ/KJzLrL/cUzGY3BAwA35L0HDANfA64DPhSzWBm/eCiYDa7EeAl4DCwlbS04QPA6oh4qmIus75wUTCb3TDwYER8rHYQsxL8noLZ7EaA3bVDmJXiomA2u/fiomBnEM99ZGZmDfcUzMys4aJgZmYNFwUzM2u4KJiZWcNFwczMGi4KZmbWcFEwM7PGfwFtLJBSTGT0JwAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "%matplotlib inline\n", "\n", "# Define the characteristic equation function\n", "def characteristic(Pe, beta):\n", " return beta * np.cos(beta) / np.sin(beta) - beta**2/Pe + Pe/4\n", "\n", "# Compute the Peclet number\n", "Pe = U * L / D\n", "\n", "# Make a list of the first few singularities\n", "singularities = [np.pi * i for i in range(11)]\n", "\n", "# Make a customized plot area\n", "fig, ax = plt.subplots()\n", "ax.set_ylim(-10,10)\n", "ax.set_xlim(0, np.pi * 10)\n", "ax.axhline(y=0, c = 'k', lw = 1)\n", "ax.set_xlabel(r'$\\beta$', weight = 'bold', fontsize = 14)\n", "ax.set_ylabel(r'$F(Pe, \\beta)=\\beta$ $\\cot$ $\\beta - \\frac{\\beta}{Pe} + \\frac{Pe}{4}$', weight = 'bold', fontsize = 14)\n", "ax.set_yticklabels([])\n", "ax.set_yticks([])\n", "ax.set_xticks(singularities)\n", "ax.set_xticklabels(['{}$\\pi$'.format(i) for i in range(len(singularities))])\n", "ax.set_title('Characteristic Equation for Eigenvalues', weight = 'bold', fontsize = 14)\n", "\n", "# Go through each interval (n * pi through (n+1) * pi) to plot \n", "# the function and singularities\n", "for i in range(len(singularities)-1):\n", " s1 = singularities[i]\n", " s2 = singularities[i+1]\n", " xs = np.arange(s1 + np.pi/100, s2, np.pi/100)\n", " ys = characteristic(Pe, xs)\n", " ax.plot(xs,ys, c = 'r', lw = 1.5)\n", " ax.axvline(x=s2, c = 'k', ls = '--', lw = 1)\n", " \n", "# add an annotation to point out the first few betas\n", "arrowprops={'arrowstyle':'-|>','connectionstyle':'arc3,rad=-.4','fc':'k'}\n", "ax.annotate(r'$\\beta_{1}$', xy = [1.54,0], xytext = [1,3], fontsize = 12, arrowprops = arrowprops)\n", "ax.annotate(r'$\\beta_{2}$', xy = [3.88,0], xytext = [4.3,3], fontsize = 12, arrowprops = arrowprops)\n", "ax.annotate(r'$\\beta_{3}$', xy = [6.72,0], xytext = [7.5,3], fontsize = 12, arrowprops = arrowprops)\n", "\n", "# make a formatted manual legend\n", "ls = [lines.Line2D([-1],[-1], c = 'r', lw = 1.5), \n", " lines.Line2D([-1],[-1], c = 'k', ls = '--', lw = 1)]\n", "labels = ['Characteristic Equation','Singularities']\n", "leg = plt.legend(loc = 2, facecolor = 'white', framealpha = 1, handles = ls, labels = labels)\n", "leg.get_frame().set_edgecolor('k')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The graph above shows the value of the function, the singularities that every $n\\pi$, and the location of the first eigenvalue near $\\beta = 1.54$. To use the model result, the first task is to identify the eigenvalue across each interval. \n", "\n", "Note that as $\\beta$ becomes larger, it tends to approach the singularity. A very good first guess for the next value of $\\beta$ each time is to assume it is $\\pi$ larger than the last value. After computing the first ~10 or so values, it is reasonable to assume the rest are just $\\pi$ more than the last. \n", "\n", "Also note that there is no way to get all the eigenvalues, since there are an infinite number. The number required for a good approximation decreases as time goes by, due to the negative exponential term. At the beginning ~100 or so is needed.\n", "\n", "### Determining the eigenvalues\n", "Scientific Python has an optimization library with a variety of methods to determine the root of a function. The example below uses [Brent's method](https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.brentq.html#scipy.optimize.brentq), which requires a function of one variable and computes the root across an interval. The code below creates a list of intervals and then uses the method to find all the eigenvalues." ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "# Compute the Peclet number\n", "Pe = U * L / D\n", "\n", "# Define the characteristic equation function with only one argument\n", "def characteristic(beta):\n", " return beta * np.cos(beta) / np.sin(beta) - beta**2/Pe + Pe/4\n", "\n", "# Make a list of the intervals to look for each value of beta\n", "intervals = [np.pi * i for i in range(n)]\n", "\n", "# Store the eigenvalues in a list\n", "betas = []\n", "\n", "# iterate through the interval and find the beta value\n", "for i in range(len(intervals)-1):\n", " mi = intervals[i] + 10**-10\n", " ma = intervals[i+1] - 10**-10\n", " \n", " # Brent's method can find the value of the \n", " # characteristic equation within a given interval\n", " betas.append(optimize.brentq(characteristic, mi, ma))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Estimating the concentration profiles\n", "Now that the eigenvalues have been determined, it is possible to evaluate the concentration at any point in time and space. The code below determines the concentration profiles across the column at different points in time." ] }, { "cell_type": "code", "execution_count": 273, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAhYAAAEaCAYAAABNbR41AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzsnXd8lEX+x9+zm7KbzaYnm00IPUCooUi1ADZsB4p6nHoellM5Pe93WMCzYEPxkBP1UJTTs+sp56FiO1QsByhEgdBJKCGk92yy6Tu/P54nYROSEJINm4R5v17Pa/eZ+T4z39k82eezM9+ZEVJKFAqFQqFQKDyBwdsOKBQKhUKh6DkoYaFQKBQKhcJjKGGhUCgUCoXCYyhhoVAoFAqFwmMoYaFQKBQKhcJjKGGhUCgUCoXCYyhh0UMRQhwWQkghxGve9kXhXYQQc/V7QQoh+nrbn85ACDFWCLFBCFGut3ONEGKqW7un6nYP16d52WWFoseihEUXQgjhL4T4sxBikxCiRAhRIYRIEUK8IoRI8LZ/nY0QwiCEuFEI8a0QolAIUSWEOCSEeFcIMdHb/nkCvW1SCPHtKSw7D/hJP6o8XW8r/sgmR5UQ4qAQ4mUhRB8PV/cKMBkQwBZgH1DKsXaXerg+hULRAj7edkChIYQIBb4GRutJZUAKEAfcCOwA9njHu85HCOEPfARcqCdVAvuBaGAOkAP86B3vvIcQwk9KWd2RMqSUnwKfesil9pAPHABCgMHA74GLhRAJUkpHcxe0o93D9NcVUsp73NJ7hCBVKLoTqsei6/B3jomKp4EwKeVIKWUoMB1NWAAghDhTCPGl3qtRJYTYJ4S4Xwjh21LhQoi+br8c57qlNxoyaWJ3txDiQyGEUwixU693tBBis97l/D8hxGC3sl7TrzsshLhKCLFXt/ve3a4FHuKYqHgbiJRSjpBSRgJjge/d6hmu+5UvhKjWezWeFkIEutk0/HoXQtyu++QQQqwVQkQ3+Wyu1tvicGvrLLf8M/Tr6ntRdgghbmhSRv1n9pQQ4u9CiAIhRK4Q4lkhhE+9DXCOfsk57sMTTYYrrhJCJAkhqtEewGOFEF8LIbL0+suFEFuEENe519/Gsvu6XfMrIcQPQogyIUSlEGK7EGKeEEKcTLvawKdSyolSyiFo9zZALHCuXkf9PfimEGKZECIfrZcBIYRZCLFYCJGq/60LhRCfCCHG6PlT9bbX+3K3XtbDopmhkJYQQswRQvyof7blQohvhBBT3PKNQojHdT8qhBBFQohtQoglbfwMFIrTBymlOrx8AMFADSCBbYBoxXaqm20RWpev1I933ewO62mv6ed93ezmttGuEjiI1nsi0XoN8vQ6q/W0DW5lvaan1ej5ewBXU7tm2iT0siWQBZhbsU0AHLptGbAbqKuvAzDodt/qadVABVrvR3273nYr7y639FIgWS/3YT1/MtrwQX37d7nZ3+VWjnSrrwA46pb2e93mR72O+rp+1A87MNfNvgrIQOuxmglcqbfxMPALUOhme8lJlt1Xt7/OLS0HOOR2/uTJtKuVv1W93WtuaU+7pc9qcg9W6ccO4Ac9b52b/R63NjqBRGCM3s56m6P6+c1o/yv16VP18h6uT2vhHkgFjri1eZJu80c9rRbYjvY/UAmkevv7Qx3q6GqH1x1QhwQ4w+2L7fkT2H6n2x0BQvW0JW7Xj9DT6r+sX9PP+7rZzHUrrzW7L9Ee+je7pa3S7R5zSzPraa+5pV2mp/2tqV0zbYp0s/nkBO1/XbcrB/roabc1U++3+nkdkKinfainZevnARwTTT8BIW7pQ/T33+j53wG+etr9HHuAm/S0+voPoglFE5o4kMB7bv7X+/Vtk3bNdSvjbY4JJCOaOLC52ZrQRIcE3jzJsvvqaWn6+Ra9PAG8y7EHatjJtKuFv1X9tXloD/u9bmkZgLXJPVgFjHRr9zQ3+7v19Gg0QS2BfzdT18NuaVPd0qfqaQ/XpzVzDzyhpxnQ7n0JrNPTntfPX3Er3wxM9vb3hzrU0dUONRTSNRBu7+UJbM/QX7+QUhbp799xyx/nMa+0LmyJ9sVfzyf660G3tKgm15VIKevtdrdiV0972r9BSpmmv2+t/TuklNua+FLvxzDAor9/QUpZDCCldEop9+rpE/TXs4Fqvdv9cT3NyrGx/Xo+llKWSCkr0XoBAGwnaFNTnpdSunRf6tB6fZYJITKFELVoPTADdduYkywbIUQU0Fs//Y+UslL/O7+rp/kCo5pc1pF2RaB9jv3R7qVVaA/kpvEV66WUydDQ7jPc8t7R07OB9XqaJ+5193vgPv3vWwdcoKfVx2isRbs3b9SHpL4DFqOCQhWK41DBm12DfWhdrD7AmUIIoX/Rt8aJ8luzN7q9D27lmvovzdpm0tzLcxcGAMVu72tbsasnD8hFe+CPE0KY9AdYa7S1/c350pIfrZEJpDeT7uqE+rKbnL8FnMex4QAHMBRN2BjpGKfic3xdSjm3DXZN2+3Oyd7vbcW9DXuBkubqlVJ+qcd1XIUmukajic3fCyGGSimbuzcUitMS1WPRBZBSlgDv66ejgSfcA+OEEGcLIabrp1v014uENpME4Bq34pJaqCbX7f0AvdypaJH6XkUXUf/QT+3AS0KI+l+RCCFGCSEu10/r2z9FHJuy2Jb2N8cutCEVgNuEEEF6fSZxLNi0vr5M4FypBSFOBC4Dlkspt55EfaDFBsCxX8ltof5X8yop5TDgYrTu+3aVLaXMRRtKA7hCb69Am30DWozM9pPwr7PY4vb+WgA98HaannYyf+uW2Mmxz+0btJ6U+r/xXGCRXu9IIFdKeb+U8lKO/U0CgfEe8EOh6DEoYdF1uAOof0gtBAqFEMl6hPx3wEg9bxHaL8Y44KAQYh+wQM97T0q5g2aQUlYAm/TTu4UQ69GGNZr+4vYWjwL/1d9fD+Tp7c9GC2g9R89bgvZQtQC7hBC7gBV63kZOYlqllNKJ/uBAe1AcFUJsRxNhv9HTH0B70I4DsoQQW4UQR9B+XbdnRkD9EMs4vX1ftOGaZP31Zr29B9BiHTpS9v31tmjDEwc51uZlUsrCNvjVqUgp1wNf6adLhRB70Hr3QtACJx/zQB1O4BH99A9Ahv43zkXrHaoXrVej3R9HhBA/owXRgjZssqujfigUPQklLLoIerzEZLQI9c168iC0bu/X0R+6Uspv0X6x/Rft79cPbcbDg2gP5NaYC/yAJkxi0cRMl+jClVJWARehBYp+hxbINwTtob4a+JdutweYBPxHtxmE1oZlwIX1sQknUe8y4NdoosTgVt52Pf9/wFloY+y1aEMQoAmYB9vR1KfRHpZlwAjaFicwFy2uoBIt2PD/OCY22lW2lPIttBknG9CGVOx6mX8A/tKmlpwafgU8gSZ8BqAJ4bXAFLfYmQ4hpfwrWo/Ij0AQ2j1QjPZ/V9+T9h3wGdrQyXC0YcuNwGy3eByFQoE+rVGhUCgUCoXCE6geC4VCoVAoFB5DCQuFQqFQKBQeQwkLhUKhUCgUHkMJC4VCoVAoFB6jxyyQFRERIfv27ettNxQKhaJb8fPPP+dLbbM/hcIj9Bhh0bdvX5KSPLFejkKhUJw+CCHSTmylULQdNRSiUCgUCoXCYyhhoVAoFAqFwmMoYaFQKBQKhcJjKGGhUCgUCoXCYyhhoVAoFAqFwmOcEmEhhHhVCJErhNjplhYmhFgnhEjRX0P1dCGEeE4Ikarv0DjmVPioUCgUCoWi45yqHovXgBlN0hYCX0sp44Gv9XPQdriM149bgBdPkY8KhUKhUCg6yClZx0JK+b0Qom+T5JnAVP3968C3wAI9/Q2pbbv6oxAiRAhhl1JmdYZvVYdLqNxfBEIA+ov+HlF/uOUhjqXXv3e/rmmeANHcNbq9cK/H/ZpGPjR5f1z+sWsb+3HM5rjy6k8Nx967Xyea2gq38pu2of76hrKa+NKobDdfFAqFQtHj8OYCWbZ6sSClzBJCROnpsUC6m91RPe04YSGEuAWtV4PevXu3y4nqNAeO9emgdo8/tbiLNl2MHCfO6kVLE4Eiml7brNhqRtwYjtkfL5z065uUe7zQPEH9TUSaMLQgChv55e5Tc+1vYtdI3DXnv1tdHJ/WqA1NxWxLwrG5tjXzGWjZTdrXVIC2mudWN8372eg+qP+7Nim30WdEkzLq3yqBq1B0Cl1x5c3m/tubfexLKV8GXgYYN25cu6SB9ZxeWM/pVV/esZoavdfzGjzR86RuV//ibuNmJ9t8zbH3DXnN+NPgZ5NrZZN6G9Ldy2/JhybXyWZt3dLqbV3u5TcuSzZNq/eR431p+IxcTT8H2aRckK5mPpeG9rt9Jq7m2t6CXX1d9fXr9q6m9bv71PRzbuYzblyXfj+43K91u0ea2Cmxe4poRXw0iBtat6l/f5xWaZLfbL3N5Itm7ZorRCNoRl8so6OazVMoTjXeFBY59UMcQgg7kKunHwXi3Ox6AZmnwqFGvzib/AOr3zYKb3CcsGpO2LiaEYlN7VwNBTaxa0HMtiRk2yKuOOZTs+XQpKwm540ELs2nNxap7u1okt7o+mN2ssk5zVwj3c8b6mnhuvr2ueN+2lreicpp6Ro3jEF+LWcqFKcYbwqLj4HfAUv014/c0u8QQrwHTABKOiu+QqHo6jQanqhP85o3CoVCcWJOibAQQryLFqgZIYQ4CixCExTvCyFuAo4AV+nmnwEXA6mAE7jhVPioUCgUCoWi45yqWSG/aSHr3GZsJXB753qkUCgUCoWiM1ArbyoUCoVCofAYSlgoFAqFQqHwGEpYKBQKhUKh8BhdcR2LU8q+Tf9j57frCAwNIzAsnMDQMCyh4VjDwrGEhhEQHIzBYPS2mwqFQqFQdAtOe2Hhqq2horSUvLRDOIuLkdLVKF8YDFhCw7CGhmvCQz+sYeEEhkc0nPv6+XupBQqFQqFQdB1Ei4uxdDPGjRsnk5KSOlSGq66O8pIiyouKKCss0I6iQsoK83HUnxcWUF3hPO5akzUIa1g4Vl1sWMMjm7wPx9ff1CH/FAqFwtMIIX6WUo7zth+KnsNp32PhjsFoxBoWgTUsAgbEt2hXXVlBWWEBjoL8BrHhKMjHUZhPWUEBWan7qSgtOe46kzUIa3iEfmjCIygikqCIKKwRkQSGhmEwqmEXhUKhUHRflLBoB34mM2ExvQiL6dWiTW11tS448jTRUZB/7H1eLpl7d1NZXtboGmEwNIiOoMgogiIiG7+PiMTPZO7s5ikUCoVC0W6UsOgkfPz8CIm2ExJtb9GmurICR34epfl5+msupfl5lOblkrF3F3sL8pGuxjEfpkBrQw9HUGQkwZE2TXhE2giOsmGyBHZ20xQKhUKhaBElLLyIn8lMeK/ehPdqfst3l6uOssJCTXQUaIKjXoAUZ2dyZOd2aiorGl3jH2A5JjT016CoqAYBooSHQqFQKDoTJSy6MAaDUY/BiCS2mXwpJZXlZZTm5lCal0tJbjYlebmU5uVQkpPFkR3bqKmqbHSNu/AIsdkIjopuOIKiotTsFoVCoVB0CCUsujFCCMyBVsyBVmz9Bx6XL6WkssyhiY68HEpzcxqER3F2JmnJW6mtrmp0jSU0TBcamugIsUUTHGkj2BZNYGgYwqDWVFMoFApFyyhh0YMRQmC2BmG2BrUoPJwlxVpPR042Jbk5FOdmU5KbzdE9O9nzv2/BbTqy0cdHG2KxaYIjxGYn2GbXxIctWvV2KBQKhUIJi9MZIQSWkFAsIaHEDEo4Lr+2pgZHfi4luTmU5GZTnJNNaW4OxTnZZO7bc9x6HoGhYbrQ0MVGtL1BgJgCrQghTlXTFAqFQuEllLBQtIiPry+h9lhC7cdHeEgpqXCUNgiOkuwsinOyKc7JIi35F3YVFTay9w+w6D0d2kyZ0OgYQuwxhEbHEBAcokSHQqFQ9BCUsFC0CyEEAUHBBAQFYx84+Lj8mqpKbWglJ5uSnCyKczThkZd2kNQtm3DV1TXY+prMDWIj1B6jiQ8lOhQKhaJbooSFolPw9TcREdeHiLg+x+W56uoozculKDuT4uxM7TUrs1nR4Wc2E2KL0YSHPYaQ6JgGEaJEh0KhUHQ9lLBQnHIMRqPb4mFjG+XV1dbq63RkUZSlCY/i7ExyDx8gZfPGRguGNYgOvXejYYgl2q5Eh0KhUHgJJSwUXQqjj482JBIdQ7/EFkRHViZF2VkNvR25h1JJ+WnD8aIjWuvhaBAd9lhC7TGYrUFKdCgUCkUnoYSFotvQSHQ0yaurrdXX58iiKDtT6+3IySL34PGiw99iaQhKDbXHNHqv9mJRKBSKjqGEhaJHYPTxaRAIzYkOLZA0k6LMTIqyMijKytDW6vhhfSNbS2hYY7ERrb0PtkXj4+t76hqkUCgU3RQlLBQ9HqOPD2ExsYTFxMLoxnk1VZUU52RrYiMzg6IsbXgldcuPVJSWNNgJYSAoKspNbBwTH9aICAwGtd29QqFQgBIWlOZXUFpQiX1gMEajWq76dMPX30Rk775E9u57XF5lWRlF2brYyDrW05Gxd3ejzd+Mvr6E2OzHDauE2mNVEKlCoTjtaJOwEEIYgHggFCgGUqSUda1f1T3YuymLLZ8exs9kJG5oOH1HhtNnWDhmq5+3XVN4GVNgIPaBg49bp0NKSXlxkS40MhteCzMzOPhLEq662gZbP7OZUHusFkTaIDhiCI2OxRSodppVKBQ9DyHd9oI4LlOIc4FbgRmAxS2rHPgSeFFK+U2nethGxo0bJ5OSkk76uurKWo7uLSJtRz6HdxbgLKkGAdH9gugzIoJ+IyMIi7GoX52KNuFy1eHIz6MoM4NCt16O4uxMSvJyG+29YrIG6SLjmOgI0YdZVBCp4lQhhPhZSjnO234oeg4tCgshxDfAOYAAKoBUoBQIAgYCZkAC30kpp58Sb1uhvcLCHemS5KU7OLyjgLQd+eSmOQAIijDRd2QEfUdGEBMfooZMFO2itqaGkpxsfdZKhjZtNiuDouxMygoLGtk2BJHWT5mN0WI7Qmx2fPxUb5rCc7RTWEwGHgOCO8ElRdemBHgQ2NiSQWvCohR4DXgX2Ow+9CGEMALjgWuA30kpgzznc/vwhLBoSnlxFYd35HM4OZ/0vUXU1bjwM/vQZ1gY/UZF0nt4OP7m0z5MReEBaiorj61EWh/ToQsQ9yBShCAoIvLY0Ep0DKExmgAJirRh9FH3o+LkaKew+BolKk5nSoBzW8psTViESCmLT1R6W+06m84QFu7UVNWRvqeQw8n5HN6RT4WjBoNR0GtwKP0SI+k3MgJLiNo2XOF5qpzljWI56odWirIyqXKWN9gZjEaCo2yNREdItLbbrDUiUokORbO0U1h03petorvQ4j3TaoxFawghrgLsUsrn2uuVJ+lsYeGOyyXJOVjCwe35HNqWR0meNkMgqm8Q/RMj6J8YSWi05QSlKBQdo36H2aZio354pbaqqsFWGAwERUYRHBXdsJV9iM2u7zgbjZ85wIstUXgTJSwU7aRThMUmYLyUsktM4D+VwsIdKSVFWU4Obs/j0La8hriMULuF/okRDBgdRURcoAr+VJxSpJSUFRVQkp3dsLNscU6WttNsbg6VjtJG9uag4AbBEWyzu72PxhISqu7fHowSFop20nWFhRDiz8DNaIGgO4AbADvwHhAG/AL8VkpZ3Vo53hIWTXEUVnJoex4Ht+WRub8YKcEaZqL/6EgGjI4kun8wwqC+pBXepbK8jJKc7MaCIyebktxsHPn5SHlsCXQff39CoqKPExwhtmiCIqMw+qgVSbszSlh4n/Hjxw+eM2dOwfz58/O97ctJ0OI949VBVyFELHAnMFRKWSGEeB+YA1wMPCOlfE8IsRK4CXjRi662GWuYiZHT4hg5LY6KsmoOJ+dzYGseO747yvav0wkI9qN/oiYyYuJDMKgZJgovYLIEYuo/EFv/gcfl1dXWUJKbq4uNY4KjODuTtOSt1Fa7DbEIA9aIyOMER3CUJjpMgVbV26E4pcTGxo5YsWLF4VmzZjk8XfbKlSvDHn300diioiKfKVOmlL799tuHbTZbs2s6bdy40fz73/++78GDB039+/evXLVq1eHJkydXNGfb02hVWAgh7mwlO9qDPpiFEDVAAJAFTEebcQLwOvAw3URYuGMO9CNhcgwJk2Oorqjl8A5NZOzdmMXO7zIwBfrSf3QkA8dGEatEhqKLYPTxPbYEehOklJQXFR4THA2v2aRs3khFkyEWHz9/rBGRWMMjCNJftfPIhnNff9OpappC0W6SkpJMd911V5/Vq1enTJ482Xndddf1uemmm/qsXbv2YFPbyspKceWVVw689dZbc+699968ZcuWRV555ZUDDx48uNNkMrVvmMAD1NTU4HsK9jw60QJZLrQhimazAemBoZA/AYvR1sr4L/An4Ecp5UA9Pw74XEo5vJlrbwFuAejdu/fYtLS0jrhyyqipquPIrgIO/JLLoR0F1FbVYbb60n90FPFjo7DHh2BQwyWKbkiV06kNrehDKo6CPBz5eTgK8iktyKO8uKjRImGgLRTWSHiER2KNiCQoPBJrRASBoeEYjF0ilKtH0pOGQmbNmtXv448/DvPz85MGg0HOnz8/8/HHH8/xRNl33HFHbFpamt8nn3xyCGDXrl3+iYmJw3Jzc7eFhoa63G0//PDDoNtuu61vdnZ2ssGg/WC02+0jnn/++bQrr7yytGnZ48ePHzxp0iTHDz/8ELRv3z5zYmJi2erVqw/Z7fZagLfffjt40aJFvXJycnwTEhIqVq5cmTZmzJhKACHE2B07duwcPnx4FcDs2bP7xsbGVj/33HOZa9eutd500039br755tyXXnrJduaZZ5a++OKL6ddcc03fpKSkQIPBwMCBAys2b968z3jy/2PtHgr5npaFRYcRQoQCM4F+aEuFfwBc1Ixpsz5IKV8GXgYtxqKT3PQ4vv5GBoyJYsCYKGqq6ziys4DUn3PZ92MWu77PwBzkx8DRkQwcF4V9QIiKyVB0G/wDArD1G4Ct34Bm8+tqaygrLMCRrwkNTXRowqMkN4eje3ZSVV7e6BohDFjCwjSh0UyPhzUiErM1SA25eIl7Vm+P25/t6NRpRYOirc6lV45KP5HdmjVrDsXGxga2NhSSkpLiN2bMmKEtlbF06dIjt912W2HT9D179pgmTpzYcHMOGzasytfXV+7cudN01llnOd1td+zYYRoyZEhFvagAGDJkSMWOHTvMzQkLgA8//DDs008/Tenfv3/19OnTBz322GO2F154ISM5Odn/5ptv7v/OO+8cuPjiix2PPfZY1KxZswbu379/V1t6PwoKCnwLCwuN6enpyXV1dSxYsCDGbrdX5+fnbwdYv369xdP/OycSFhecKGiyg5wHHJJS5gEIIT5EW9EtRAjhI6WsBXoBmZ3og1fx9XMTGVV1pO0sIPXnHPZszGLHdxkEBPsxcEwUA8dGqcBPRbfH6ONLcJQWg9ES1RVOrYejXnTU93jk55FzMJXULZuoq61tdI2Pn7/e2xFBYGgYASGhWIJDsISENrwPCAnFHGhFGNSQ4+lMfHx8tcPh2Hay1zmdTmNwcHCjeIrAwMC6kpKS437ql5WVGYOCghrZWq3WOofD0WK3wG9+85uCkSNHVgFcccUVhZ9++mkIwJtvvhk2bdq0kssvv7wU4JFHHsl5+eWXbV999VXgpZdeesI4EiGEXLZsWabZbJYAvr6+MicnxzclJcVv+PDhVTNmzChrS/tPhhMJi0IhxDrgU+BTKWWWh+s/AkwUQgSgDYWci9bFth64Em1myO+Ajzxcb5fE19/IwLGaiKiurCVth9aTseuHTJLXHyUw1J8Busiw9VO/0BQ9Ez9zAOG9ehPeq3ez+dLlwllagqMgv6HHo1QXH478PI7u3UV5cRF1NTXHXWswGgkICiYgOBRLSMhxr+agYMzWIEyBgZgCrfj6m9T/2QloS09CTyAgIKCutLS0kSotLy8/TmyAJjgcDkcj27KyMqPVam1x887o6OiGGzYgIMDldDoNAJmZmb5xcXENP/CNRiN2u706PT29TcESoaGhtQEBAQ09G4sWLcq+9957Y2bMmDEI4Prrr8974oknsttSVls5kbC4Hm2GxiPAS0KI7cBaNJHxU0crl1L+JIRYjTaltBbYija08SnwnhDicT3tlY7W1d3wM/kQf4aN+DNsVFfUcig5n9SfcxtmlwSG+TNwrI1BZ9jUOhmK0wphMGAJCcUSEkr0gPhmbaSUVFc4KS8uwllcTHlJEeXFxTibvOalp+EsLm60I607BqNPg8gwBVoxW62YLNZGaabAQMyBQfhbLPiaTPiZzPiaTPj6m9Rqp12QlJQUv1GjRg1rKX/ZsmVp8+bNO24oJCEhoTI5OblhyGf37t1+1dXVYvjw4ZVNbUeMGFH5wgsvRLtcLuqHQ/bu3Wv+wx/+kHuy/sbExNTs2rWrYVdCl8tFVlaWX1xcXA2AyWRylZeXN4iY3Nxc39jY2AYh0vTZEBoa6lq1atVR4GhSUpLpggsuGDxhwoTymTNnemwWTat3vZTyQ+BD3bmxaCLjEuAvQogC4HM0ofGllLJdTkkpFwGLmiQfRNuLRAH4mX0YPCGawROiqaqo5fD2PFJ+ziX563S2rTtCaHRAgwgJiVIrKCoUQgj8Ayz4B1gIi+nVqq2UksryMpzFRVSUllJR7qCyzEFlWZn26tDPyx2U5ueRe/gglQ4HNVXHPU+Ow+jjg6+b0PAzmZo/10WIwccHY8Phi8Fo1NN9Mfr4IAwGhEFgEAYQQjuA8Ng4rOERHvnsujsRERE1qamp/kCzz6T4+Phqp9O59WTLnTt3bsHUqVMTvvjii8DJkyc777vvvtgLL7ywuGngJsDFF1/sMBqNcvHixVF333133jPPPBMB0Jahi6Zcd911hRMnThz60UcfWWfMmFG2ePHiKD8/P3neeeeVASQkJFS8/vrrYWPHjs1Ys2ZN0JYtW6yjR48ub6m8d999N3jEiBGVQ4cOrQoJCakzGo2yHYGbrdJmOS2l/Bn4GXhMCBGFFmR5MbAKeAZ41KOeKZrF3+zD4Il2Bk+0U1lWQ+ovuaRsyWHzJ4fY/MkhovoGMegMGwPHRWEJVnt9ESCxAAAgAElEQVSXKBQnQgiBOdCKOdB6UtfV1tRQVa6Jj4oyB1XlZdRUVlJdWUlNZSU1lRXUVDU+r67SXssKCqipqtDzKqiprGq0KNnJct7Nf2DU+Re3+/qexD333JN9zz33xD3yyCO9/vznP2c9+uijHpkVMm7cuMqnn3467YYbbuhXXFzsM3ny5NJ33nnncH3+2WefHT958mTHkiVLsk0mk/zggw9Sb7nllr6LFy/u1b9//4oPPvggtT1TTUeNGlX10ksvHZo/f37v6667znfIkCEVa9asSakva/ny5UduvPHGfsHBwVHnn39+8fnnn1/UWnn79+/3v/vuu3sXFhb6BAUF1c2dOzevPYKnNdq98mZDAdpOp2H1AZjeoqusvOktHIWVpCblsn9LNvnpZQgBsYNDGTTeRv/RUWoXVoWii+Ny1VFXW4urto662hpctbXU6Yerrpa6mhqklEiXC5fL1TBtVyIJsdkJDA1rV709abqp4pTSsSW99dkazVEF7AP+KaX06iISp7uwcKcwq5yULTns35JDaV4FRh8DfUaEM+gMG31GhOPjq9YEUCgUGkpYKNpJh4VF/UJZ7lEg9ecSKAPOkVKe9BQeT6GExfFIKck97GD/lmxSknKpKK3Gz2Sk/+hIBp0RTeyQULUQl0JxmqOEhaKddHivkL8D89Bma+wGEoDLgH8AfYHz0WIsftURLxWeRQiBrV8Qtn5BTJk9kIz9xezfksPBX3LZuykbS7Afg8ZHM3hiNOGxgd52V6FQKBQ9gLb2WKwFnFLKq93S3gfMUsrL9KGSCVLK4zcXOEWoHou2U1tTx+HkAvb9lM2RnQW4XJKIuEAGT4gm/gybCvpUKE4jVI+Fop10uMdiGrBLCGGSUlYKIfyB3kD9/h1baX4pbkUXxMf32EJcFY5qUpJy2PdjNhtWp7Lx36nEDQ1j8IRo+iVG4uun4jEUCoVC0XbaKix2AGcAOUKIdLRltq1A/SJZ44CjnndP0dmYrX4N27wXZZez78ds9m3OZt2ru/U9TSIZPCGa2EGhajlxhUKhUJyQtgqLG4BPgP5A/eYtqcCNQogg4AD6QlqK7ktotIWJswYw4Vf9yUwpZt9P2aTq8RiBof5aPMaEaMJiLN52VaFQKBRdlDavY6GvVzERiAUy0LY2b3Hd81NNe2MsKnbuouTjj4iYNw+f0NBO8Kx7U1Ndx+Ht+Vo8xu5CpEsS2dvKkEl2Bp1hwxTYpuXqFQpFF0XFWCjaScemmzYYCxEONPq5KqU80n6/PEd7hUXhG2+Qs+QpDBYL4bf8nrDf/haDydQJHnZ/nKXVpGzJYe+PWeSnl2EwCvqNjGDIJDu9h4VhMKpdIxWK7oYSFt5n9uzZfWNjY6ufe+657rSTd4v3TJueBEKI84UQR4Bc4JDbcdAj7nmRsOuvp/9HawgYO5a8ZX/jwEUXU/LRR0hX+5fX7akEBPkx6tw4fn3/eH79wBmMmNqLzNRiPn0hmdfu28iG1SkUZHh8B16FQqFoM7GxsSPWrFlzcuuzt5GVK1eGxcTEjDCbzaPPO++8ATk5Oa1tg96nb9++ww0Gw9jnnnsuvDP86aq09SfmSrSATdHk6BE/Uf3j44lb+SK9X38dn/BwMhcs5NCVV1K+aZO3XeuyRPSycuZV8fxuyRQunjcCe/9gkr85ynuPbeb9J7aQvP4olWXHb1utUCgU3ZGkpCTTXXfd1eeVV145lJ2dvd1sNrtuuummPi3Zjxw50vnss8+mDR061Hkq/WyNmppT853cVmEQDnwJBEkpDe5HJ/p2yrFMGE/f9/9FzNKl1BUXc+SGG0m/9TaqDhzwtmtdFqPRQL9RkVx02wjmPjWFM6+OR0rJD//azz8X/I8vXtrB4eR8XHWqB0ihUHQus2bN6peVleU3Z86c+ICAgNEPPPCAzVNlv/baa+HTp08vvuiii8qCg4NdS5Ysyfzyyy9DioqKmn0O3nfffXkzZ850+Pv7t+nLr6ioyGfq1KkDLRbL6JEjRw7ZtWtXw4JC69atswwfPjzBarUmDh8+PGHdunUNIQlNe2jmz58fM3PmzH4A+/bt8xNCjH3mmWci7Hb7iEmTJg12Op1i5syZ/UJCQhLry0tPT/foZlJtLex54AogVgixX3Z057IujDAYCL7sUqwXnE/Rm2+Sv/IlDv5qJqG/vpqIO+7AJ6x9G/2cDpitfoyaHseo6XHkH3Wwd1M2+zdnc2BrHuYgPwaPtzFkkl2t8qlQ9CTW3B5H7u6ATq0jaqiTWSvST+jKmjWHYmNjA1esWHF41qxZze7YmZKS4jdmzJihzeUBLF269Mhtt91W2DR9z549pokTJzZsRz5s2LAqX19fuXPnTtNZZ53V4V6Jjz/+OGzNmjX7zzzzTOfs2bP7LViwIHbt2rUHc3JyjLNnz45/8sknj9xyyy2Fr776atjs2bPj9+/fvyM6OrpNEyi+//77wH379u0yGAxyxYoV4Q6Hw5ienp5sNptdmzZtCrBYLB795ddWYfFv4A605bwRomE9Ayml7JHbZhr8/Qm/+WaCr7iC/L//naJ/vU/Jx58QfusthF1/PQZ/tTpla2hDJVYmXTGAIzsL2Lspm+RvjrLtq3Si+lhJmBJD/Bk2teuqQqE4pcTHx1c7HI6T3tfK6XQag4ODGz3IAwMD60pKSjyyiuCMGTOKpk2b5gS49tprCxcuXNgLYPXq1cF9+vSpuv322wsBbr311sIXX3wx6v333w+58847C9pS9uLFizODgoJcAL6+vrKoqMhn9+7d/hMmTKjwhChqSlu/1d8GgptJ7/ErJvmEhRH90EOEXnstuUufJm/Z3yh+9z2i7r4L60UXuYssRTPUD5X0GxVJhaOa/Ztz2LMxk+/e2ceGD1IYMCaKhCl2YuJD1GepUHRH2tCT0BMICAioKy0tbTTsUV5efpzYaC82m60hAMJisbicTqcRIDMz069Xr15V7ra9evWqzsjIaPNc/wEDBjSUPW/evML09HS/a665pr/D4TBeccUVhc8++2yGv7+/x0Yi2hoj0RvYgra09+gmx2mB/4ABWoDna//EEBRExvy7SLvmWiqSk73tWrfBbNVnlTwwnisXjmPwJDuHtuex5m9beeuhH0n6/DBlRVUnLkihUCjaSUpKil9AQMDolo4XX3yx2fHuhISEyuTk5IYhn927d/tVV1eL4cOHV3amvzExMdVHjx5t1EWekZHhFxsbWwNgNptd5eXlDc/y7Ozs4zoMDAZDg2jw9/eXy5Ytyzpw4MCuH374Ye+6deuCX3jhBY/OWmmrsHgZqAY2SSm3ux+edKY7YJk4kX7/Xo398ceoTk/n8NW/JuOee6nJyvK2a90GIQS2vkFMvWYwc/96JufNTSAwxJ+fPjrIG3/ZwNoV2zm4NY+6WhXwqVAoTo6IiIia1NTUFseq4+Pjq51O59aWjnnz5h0XXwEwd+7cgm+++Sbkiy++CCwtLTXcd999sRdeeGFxaGhos19UlZWVwul0CimlqKmpEU6nU9TVnXznxuzZs0sOHz7sv3LlyrCamhpWrVoVmpqaarrqqqtKAIYOHep87733wqqqqsT3338f8Pnnn7e60uMnn3xi3bx5s7m2tpaQkJA6Hx8faTQaPRo32VZhcS4wCcgTQmwXQvyiHz970pnugjAaCbnySgZ88QXht96K48svOXDRxeQ99xwuZ5eZWdQt8PUzMniincvvGsO1j05kzIV9yD/i4POXdvD6fRvYsDqFwqzyExekUCgUwD333JO9bNkyu9VqTXzooYc8Nitk3LhxlU8//XTaDTfc0M9ms40qKyszvPLKK2n1+WeffXb8woULo93OB1ksljFbt2613H333X0sFsuYzz///KTX14iOjq5bvXp16vPPP28LCwtLXL58efTq1atT7XZ7LcCSJUsy0tLS/ENDQxMfeuihmJkzZzYrjOrJzMz0vfrqqwdYrdbRQ4cOHT5p0iTHvHnz2hSr0Vbaum16Sz8dpZSyS2x/6c1t02syMshd9jdKP/sMH5uNqLvmE3TppQhDj5qNe8pw1bk4sruQPRuzOLw9H5dLYusXxNApMQwcF4WfSQV8KhSeQq28qWgnHVvSWwjxu5bypJSvt9Mpj+JNYVGP85et5DzxBJU7d2IaNZLo++7DnJjoVZ+6O87Savb9lM2eDZkUZTvx8de2fB862U70gGAV8KlQdBAlLBTtxDN7hXRluoKwAJAuFyUffUze3/5GbV4eQZddRtRd8/GNjj7xxYoWkVKSc6iUPRsySUnKpaaqjhBbAAmT7QyeGI0lWE3/VSjagxIWinZy8sJCCHED8KaUsrbFi4XwAa6XUr7aYRc7SFcRFvW4ysvJf3kVhf/8JxgMhP/+ZsJvukltcOYBqitrOfBLLns2ZpGVWoIwCPoMD2foFDt9hoerzdAUipNACQtFO2mXsHABeWiLY30P7AEcgBVIAM5BW40zoivEWXQ1YVFP9dEMcpcuxfHll/jE2LHdey/WCy9UXfgeoii7nD0bs9j7YzYVpdUEBPkxZFI0CZNjCLF17mKACkVPQAkLRTtpl7C4EXgYbfOx5owEcBRYJKX8Z8d97BhdVVjUU755MzlPPEnV3r0EjBuH7f6/YEpI8LZbPYa6OhdHdhawe0MWaTsLkC6JfWAwCZNjGDg2Cl9/r2tfhaJLooSFop20L8ZCCGEELgMuAUYCoUAxkAysBT6RUnpk1bGO0tWFBYCsq6P4g9XkLV9OXWkpIVddReSf7lT7j3iY8pIq9v2YzZ6NWRTnOPH1NxI/LoqEKTHY+gWp3iKFwg0lLBTtRAVvdiXqSkrIW7GCorffwWCxEPnHPxL6mzkIHzWN0pNIKck6UMKeDZmk/pxLbbWLsBhLQ8CnOdDP2y4qFF5HCQtFO1HCoitSlZpKzhNPUL5xE/7x8djuvx/LxAnedqtHUl1ZS2pSLrs3ZJJzqBSDUdBvVCRDz7QTNyQMYVC9GIrTEyUsFO1ECYuuipQSx1dfkbvkKWoyMrDOmIHt3nvwjYnxtms9loKMMvZsyGLvT1lUldcSGOZPwuQYEibbsYapWTuK0wslLLzP+PHjB8+ZM6dg/vz5+d725SRo8Z5R8/K8jBCCoPPPp/+na4m484+UffstBy6+hLwVK3BVdureNqct4bGBnHl1PDcsOZMLbh5GaLSFLZ8e4o37N/Lxc9tI/TmXuhq1T4lC0R2JjY0dsWbNmpNeOrstrFy5MiwmJmaE2Wwefd555w3IyclpMSpcCDHWbDY3bG7261//uk9n+NQV8bqwEEKECCFWCyH2CiH2CCEmCSHChBDrhBAp+murm6r0BAwmE5F/+AMDPvuUwGlTyX/+7xy8+BJK//tfekqvUlfD6GsgfpyNX92ZyG8fn8QZF/elKKucL1ft5LWFG/jf+ykUZJR5202FQtEFSEpKMt111119XnnllUPZ2dnbzWaz66abbmpVLGzZsmV3/eZm//rXv9Jasz0V1NTUnNjIA3hdWADPAl9IKYcAo9DWy1gIfC2ljAe+1s9PC3xjYuj1zDP0fv11DIGBZNz5J9Jvuomq1FRvu9ajCQo3M/6y/vx28WQu++MoYgeHsOO7o7z32GZWP5XE7v9lUl3Z4lpxCoWiCzBr1qx+WVlZfnPmzIkPCAgY/cADD3hsE7LXXnstfPr06cUXXXRRWXBwsGvJkiWZX375ZUhRUZFHnqNpaWl+Y8aMGWKxWEZPmTIlPisrqyGa/+233w4eOHDgMKvVmjh+/PjBv/zyS8OYrRBi7M6dOxuWHp49e3bfO++8MwZg7dq1VpvNNvL++++PjoiIGHXVVVf1y8rK8pk2bdpAq9WaGBwcnDh27NjB7dl1tTXaNA1BCDESeB4YDVjcsqSUst1TGYQQQcDZwFy9sGqgWggxE5iqm70OfAssaG893RHLhPH0+/DfFL33L/Kee46DM2cRdt21RNx+O8agIG+712MxGAS9h4XTe1g4FQ5tn5LdG7JY/9Zefvgghfix2rTV6P5q2qpCAfDghgfjUotSO3U1uoGhA52PTXks/UR2a9asORQbGxu4YsWKw7NmzXI0Z5OSkuI3ZsyYoS2VsXTp0iO33XbbcTuE7tmzxzRx4sSGrZaHDRtW5evrK3fu3Gk666yzmt3Wevr06YNdLpcYM2ZM2fPPP58+ePDg6pbq/fDDD8M+/fTTlP79+1dPnz590GOPPWZ74YUXMpKTk/1vvvnm/u+8886Biy++2PHYY49FzZo1a+D+/ft3mUymE3ZnFxQU+BYWFhrT09OT6+rqWLBgQYzdbq/Oz8/fDrB+/XqLp7/L2ioK3gRGNJPeUW/6o63u+U8hxCjgZ+BPgE1KmQUgpcwSQkQ1d7EQ4hbgFoDevXt30JWuh/DxIey6awm65GLylj9L4RtvUvLJWqLm/5ngK65Qu6d2MmarH4nn9WbUuXHkHCplt75PyZ6NWYRGB5AwJYbBE6IJCFLTVhWK7kJ8fHy1w+HYdrLXOZ1OY3BwcKOf9oGBgXUlJSXNxll89tln+6ZNm1ZeVlZmmD9/fuyll14av3v37l2+vr7Nlv+b3/ymYOTIkVUAV1xxReGnn34aAvDmm2+GTZs2reTyyy8vBXjkkUdyXn75ZdtXX30VeOmllzYrntwRQshly5Zlms1mCeDr6ytzcnJ8U1JS/IYPH141Y8YMj4/3tlVYDAR2AX9EWyDLU4P+PsAY4I9Syp+EEM9yEsMeUsqXgZdBmxXiIZ+6HD6hodgfeZiQq68iZ/ETZD3wIEXv/YvoB+5Xu6eeAoQQRPcPJrp/MGdeFU/qz7ns2ZDJxn+n8uN/DtBvVAQJU2KIGxqGQU1bVZxmtKUnoScQEBBQV1pa2ujXXHl5+XFio56LLrqoDMBkMtW9+uqrR6xW6+itW7eax48fX9GcfXR0dEMAREBAgMvpdBoAMjMzfePi4hp6OoxGI3a7vTo9Pb15hdKE0NDQ2oCAgIbn46JFi7LvvffemBkzZgwCuP766/OeeOKJ7LaU1VbaKiy+BmqklN96snK0JcGPSil/0s9XowmLHCGEXe+tsAO5Hq63ASklJVUlhJhCOqsKj2EeNow+b79F6dpPyV26lMNzfkPw5ZcTddd8fCIivO3eaYGfyYehU2IYOiWGwsxydm/MZN+P2RzYmkdgqD9DJtlJmGwnKMLsbVcVCkUzpKSk+I0aNWpYS/nLli1Lmzdv3nFDIQkJCZXJyckNQz67d+/2q66uFsOHD2/T9D0hRLsC8WNiYmp27drV8IXicrnIysryi4uLqwEwmUyu8vLyBsGTm5vrGxsb2yBEmg5zhIaGulatWnUUOJqUlGS64IILBk+YMKF85syZJ+z9aCtt7UvPAGYJIT4WQjwqhHio/uhI5VLKbCBdCDFYTzoX2A18DPxOT/sd8FFH6mmNt/e8zayPZrEt96R7xryCEILgyy6l/2efEX7zTZSsXcuBGRdR+MYbyFMU8avQCIuxcOaV8cxdMoULfz+cMLuFpM8P8+aDm/ho+VZSknLUtFWF4hQTERFRk5qa6t9Sfnx8fHX9TI3mjuZEBcDcuXMLvvnmm5AvvvgisLS01HDffffFXnjhhcWhoaHH/ZMnJSWZNm7caK6traWkpMRwyy23xEVFRdUkJiae9BoC1113XeH69euDP/roI2tVVZV4+OGHbX5+fvK8884rA0hISKh4/fXXw2pra1m9enXQli1bWp1q++677wbv3LnT3+VyERISUmc0GqXR6Nm9lNoqLG5Fi6e4FLgfWIS2QdkiD/jwR+BtIUQykAg8ASwBzhdCpADn6+edwqSYSQT4BnDjlzfyn5T/dFY1HscYaCHq7rvp/9FHmBMTyXniSQ5dcQXlP/504osVHsXoY2Dg2CguuzOR6xdPZvyl/SjJreC//9jFPxf+jx/+tZ+8dI/9GFAoFK1wzz33ZC9btsxutVoTH3roIY/NChk3blzl008/nXbDDTf0s9lso8rKygyvvPJKwxTSs88+O37hwoXRoA1fXHPNNQOsVuvofv36jThy5IjfJ598kuLv73/SXRajRo2qeumllw7Nnz+/d0RExKjPP/88ZM2aNSn1gZvLly8/8t///jckODh49FtvvRV+/vnnF7VW3v79+/0vvPDCQfrsk4S5c+fmtSVW42Ro08qbQojXaCGuQkp5gycdai8dWXmzpKqEu7+7mx+zfuS6hOu4a9xd+Bi6z74dUkrKvvmGnCee1FbvvGgGtnvvxddu97Zrpy3SJTm6t4jdGzI5uD0PV60kIi6QIRPtDBpvw2xVAZ+KroFaeVPRTtSS3iei1lXLsqRlvLXnLSbZJ7H0nKUE+wd70MPOx1VZScErr1Dw8iowGIi49VbCbrwBg596iHmTyrIaUpJy2Lspi9w0BwaDoM+IcIZMstNneDhGHzW7R+E9lLBQtJOOCwshxCTgdqAPcBh4QUq5yRPeeQJP7RXyn5T/8OiPjxJjieH56c/TP6S/B7w7tVQfzSD3qadwrFuHb5/e2O67D+vUqd52S4G2T8neH7PZ91M2FaXVmAJ9GTTeRsJkOxG9OmUVYoWiVZSwULSTjgkLIcSFwFrAPcLDBVwmpfy8w+55AE9uQrY1dyv/t/7/qKqr4q9n/5Wze53tkXJPNWUbNpCz+AmqDx4kcOpUbH+5D78euN5Hd8RV5+LI7kL2bsriUHK+GipReA0lLBTtpMPCYhMwFvg7sBcYAtwB/CylnOQhJzuEp3c3zS7P5s5v7mRv4V7+b+z/ccOwG7rlSouyuprCN98if8UKZE0NYTfdSMQtt2AI6NSF8hQnQatDJSPCMRrVUImi81DCQtFOOiwsSoA1UsrfuaW9DsySUnaJQITO2Da9oraCBzc8yJeHv+SS/pfw8KSHMfl0z221a3JzyX36aUo//gQfux3bgnuxXnhhtxRLPZmmQyVmqy+DzohmyORoNVSi6BSUsFC0kw4LizS0pbfPlFJWCiFMwP+ASClll9gKtjOEBWgzLv6x4x88t/U5hoUP49lpz2KzeGwG0ynH+fPPZD++mKo9ewiYOJHo+/+Cf3y8t91SNEENlShOFUpYKNpJh4XFP4AbAQeQDsQBgcA/pZQ3e8jJDtFZwqKeb458w30/3IfF18LyacsZGTmy0+rqbGRdHcXvv0/u8mdxlZVpm5vdcQdGq/pF3BVRQyWKzkQJC0U76bCwiAA+ASa4Jf+EFryZ32H3PEBnCwuAlKIU/vjNH8lz5rFo8iJ+NeBXnVpfZ1NbVETe8mcpfv99jGFhRM2fT/Dls9TmZl0YNVSi8DRKWCjaiUemmwq9oL5o002TZBdaBONUCAuA4spi7vruLjZnb+b6odfz57F/7laLaTVHxa5d5Dz2OBXbtmEaNZLoBx7EPGK4t91StIIaKlF4CiUsvM/s2bP7xsbGVj/33HOZ3vblJGjxnmnxp6kQorcQIqT+PdrwRw5aT0UOEKenn1aEmEJYef5KfjPkN7yx+w1u//p2SqpKvO1WhzAPG0afd97GvuRJajIyOXz11WQ9+CC1hc0uma/oAhiMBvqOiGDGLSO4YcmZnD1nEAaD4H8fpPDagg189mIyB7flUVer9ipRnF7ExsaOWLNmjce779LS0nynT58+MCoqaqQQYuy+fftaVe/79u3zmzBhwiCz2Ty6X79+wzrDp65Ka33eh4AH9feH9fOmx8HOdK6r4mvw5S8T/sLDkx5mc/Zmrv3sWg4Wd++PQhgMhMyaxYAvPifsd7+j+D9rtM3N3nobWVvrbfcUrWAK9GXE1F5cdd8ZzHlwPCPPjSP7UCmfr9zBPxf8j/Vv7iF9byEuV5fpYFQouh0Gg0FecMEFJe+8886Bttj/+te/7j9ixAhnXl7etkWLFmX89re/HZCZmenV7u2aU7RRZWvCQuhH03P347QejJ89aDavXPAKjmoH1352Ld8f/d7bLnUYY2AgtoUL6P/RGszDh5Hz+OMcmn0lzi1bvO2aog2ExwYyZfZA5j45mUtuH0mfYeGkJOXy8fJtvLZwA9+/u4/M1GKkEhmKHsisWbP6ZWVl+c2ZMyc+ICBg9AMPPOCxKXxxcXG1CxcuzDvnnHPKT2SbnJzsv3v37oClS5dmBgYGyrlz5xYPGjSo4q233gpt6ZqioiKfqVOnDrRYLKNHjhw5ZNeuXQ07tK5bt84yfPjwBKvVmjh8+PCEdevWWerzmvbQzJ8/P2bmzJn9QOs1EUKMfeaZZyLsdvuISZMmDXY6nWLmzJn9QkJCEuvLS09P96jgabEwKaWhufeKxoyxjeG9S97jT+v/xB1f38GdY+7kpuE3dfv1IfwHDCDulVdw/HcdOU8tIe231xN0ySVE3XsPvrbuO932dKF+qKTviAhqq+tI21lASlIOuzdmseO7DAJD/RkwNor4cTai+li7/f2q8B6Zf7k/riolpVNX3POPj3fGPLE4/UR2a9asORQbGxu4YsWKw7NmzWp2x86UlBS/MWPGDG2pjKVLlx657bbbOjQOvG3bNnOvXr2q3LdUHzZsWMWuXbtaXAjp448/DluzZs3+M8880zl79ux+CxYsiF27du3BnJwc4+zZs+OffPLJI7fcckvhq6++GjZ79uz4/fv374iOjq5riz/ff/994L59+3YZDAa5YsWKcIfDYUxPT082m82uTZs2BVgsFo+OmbZJpQghvgE+kFK+6JZ2OTBVSvknTzrUHbEH2nn9otdZtGERz/7yLPsK9/HolEcx+5i97VqHEEIQdOEFBJ59FgWr/kHBP/6BY/16Im67jbC5v1Obm3UTfPyMDBgTxYAxUVRX1nI4OZ+UpFx2rD/K9q/SCYowMXCcjfhxUYTHBiqRoejRxMfHVzscjm2dWYfD4TBYrdZGD/3g4OC6zMxM35aumTFjRtG0adOcANdee23hwoULewGsXr06uE+fPlW33357IcCtt95a+OKLL0a9//77IXfeeSnC2bMAACAASURBVGdBW/xZvHhxZlBQkAvA19dXFhUV+ezevdt/woQJFWeddZazve1sibZ2f0wFmv4hpqFtSnbaCwsAs4+Zp85+isFhg3n2l2dJK01j+bTlxATGeNu1DmMwm4m8848EXz6LnCVPkfe3v1H879XYFiwkcNpU9SDqRviZfBg0PppB46OpLK/h0PY8UpNy2frfI/zyRRqh0QEMHBtF/Bk2QqMtJy5QcdrTlp6E0w2r1eoqKytz31uL0tJSQ2BgYIs9DDabrSEAwmKxuJxOpxEgMzPTr1evXlXutr169arOyMhoUaQ0ZcCAAQ1lz5s3rzA9Pd3vmmuu6e9wOIxXXHFF4bPPPpvh7+/vsfHRVoc4hBCvCiFe1U8vqD8XQrwGzAEqPOVIT0AIwU0jbuLv5/6do46jzFk7h6TsnjMryy8ujrgVfydu1SqE0Yejf/gD6bfcStXBQ952TdEOTBZfEibHcNmdidzw1BTOuWYwAUF+bPnsMO88/BPvPb6Zn784TEme+jdX9BxSUlL8AgICRrd0vPjii2EdrSMxMbHi6NGj/kVFRQ3P2F27dgUMGzas8mTLiomJqT569Ki/e1pGRoZfbGxsDYDZbHaVl5c31JOdnX1ch4HBYGgQDf7+/nLZsmVZBw4c2PXDDz/sXbduXfALL7wQfrJ+tcaJYifmAr8DJDBUP58LXA9EAN0/WrETOLvX2bx9ydsE+wfz+//+nn/t/Ze3XfIogWedSf+P1hC1cAEVW7dy8Fe/Iuepv1LnaHZIU9ENMFv9GH52LLPmj2Huk1M48+p4fP0M/LjmIG89uIkPntzCtq+OUFZ00t+LCsUpJSIioiY1NdW/pfz4+Phqp9O5taVj3rx5LcZXOJ1OUVFRYQCorKwUTqez2e7akSNHVg0ZMsS5YMGCGKfTKd54442Qffv2ma+77rqik23P7NmzSw4fPuy/cuXKsJqaGlatWhWamppquuqqq0oAhg4d6nzvvffCqqqqxPfffx/w+eeftxggCvDJJ59YN2/ebK6trSUkJKTOx8dHGo1Gj0Zzn0hYPAI8ijYD5Cf9/BHgIeAG4EpPOtOT6Bfcj3cueYdJMZN4/KfHeXjjw1TXVXvbLY8hfH0JnzuXAV9+QfCsmRS+9hoHZlxE8b8/RLrU2gndGUuIP6OmxzH73nH8dvEkJl0xAClhw+pUXr9vIx8+/TM7vj2Ks7Tn3M+KnsM999yTvWzZMrvVak186KGHPBppbrFYxgQHB48GSExMHG6xWMbU511zzTW9r7nmmoa1nd5///2D27Zts4SFhY1+6KGHer355psHYmJiTnrufnR0dN3q1atTn3/+eVtYWFji8uXLo1evXp1qt9trAZYsWZKRlpbmHxoamvjQQw/FzJw5s9XA08zMTN+rr756gNVqHT106NDhkyZNcsybN69NsRptpa1Lei8CfpRSfunJyj3JqVp582Spc9Xx921/5x87/sGoyFE8M/UZIgMive2Wx6nYuYucx/XVO0eMIPr+v2BOTPS2WwoPUpzrJDUpl5SkHAozyxECYgeHMnBsFANGR2EKbPOQr6ILoVbeVLQTjyzp3RuYDETitr6FlPK5jnrnCbqqsKjny8Nf8uCGBwn0DeSZac8wKnKUt13yOFJKSteuJfevS6nNyyN45kwi75rP/7N33+FRVOsDx79nN5u66YU0QipFOtIUEFGKIgISUEEUL96fykVFUIpX7xULdhS5YC9gQURAmqJYsFBEVFBATMVQEkhI78nunt8fswkBEgjJht0k5/M888zu7OzMO9lk582phqAge4em2Fh2elF1kpGfWYpOJ2h7iR+xvYOI6h6Ii1vzHua+NVGJhdJAjZ6ELB54Hzir3kpKqT/7HRefoycWAAk5CczYOoPMkkz+0/8/3BB3g71DahLmomKyX3+dnGXLEAYDAf+ahu9tt6nuqS2QlJKTR4pI+uUEyb9kUphTht5JR0RnP9p18addF3+MvnV23VccgEoslAZqdGKxF+gMZAHBwJ9AHLBdSnmVjYJslOaQWIA2idmcH+awM2MnN3W4ibl95mLQt8wi5Iq0NE48+xxF336LoV0EbR56CM8rr7R3WEoTkVJy4lAByb9kkrInk6JcrYecX6gH7Tr7E9HFn5AYb/ROarw9R6ISC6WBGp1YlAHrgHRghpRSL4T4Cm2G04dsFmYjNJfEAsBkMfHyby+z7MAyegX14oXBL7TIdhdVin7cxomnn6YiNRWPKwbRZu5cXGJi7B2W0oSklORmlJC2P5u0A9lkJOdhMUsMLnrCO/rSros/EZ398fRTpRn2phILpYEanVgUAe8CucDDQCywCLhCSnnOri0XS3NKLKp8nvo583fOx8PgwcLBC+nVptf539RMycpKcj78kJNLlmIpLcX3phsJuOcenPwa3WVcaQYqykwc/SuXwwe0RKMo51RpRkRnf9p19iMk1keVZtiBSiyUBmp0YpEI/AGsB5YDJkAPnJBSOsTQks0xsQBIzE1k5taZpBel82CfB5nUcVKLHsnSlJPDySVLyf34Y3RubgTcfRe+t96KzqXObudKC1NdmnEgm8MHsklPzsNikji56AnvUFWa4YeXf/MeEr+5UImF0kCNTiweBNoBc9HGs+gMWIB7pJSv2SjIRmmuiQVAQUUBD//4MN8d/Y6RUSN59LJHcTc06Zw+dleekkLm8y9Q9N13GMLCCHpgFp7XXtuikyqldhVlJo4l5JJ2IIfD+7MpzNEG4fINdieiiz/tOvsTGuuD3qBKM5qCSiyUBmp4YiG0b/q2QLmU8oQQwoA2CudJKeUxm4bZCM05sQCwSAtv7XuLJXuWEOsby6IrFxHhFXH+NzZzxTt3cuLZ5yj/6y/cuncnaN5c3Hv2tHdYip1IKck9XqJVmeyvpTSjsx8Rnf3xClClGbaiEgulgRqdWJQDK6SUt9s2Lttp7olFle3HtjP3x7lYLBaeHvQ0g9sOtndITU6azeSvW0/WokWYsrLwvPYagh54AOfwcHuHpthZRZmJY4l5HLY2Ai3MrlGa0dlamhGnSjMaQyUW9hcfHx8ZFhZWsXjx4nR7x3IB6vydOe9fo9Qyjz2AQzTSbOkGhA3g41EfE+4Zzj3f3sOSPUswW+qcEK9FEHo9PvHjiPnyCwKmT6fou+9JvXYkJ557HnNBgb3DU+zI2dWJqG4BDJ7UgVufvIxJ8/sxYHwsRl8X9n1/lA2L9/LWAz/w2dLf2ffdUQpOqgnTWruwsLCu69at87T1cVeuXOl96aWXdvD09OwREBDQ/eabb25Xc5KxMyUkJDj369evvZubW8+oqKjOTRGTo6rv8HhbgQeFEO8CO9BKMACQUr7XFIG1ZmHGMN679j2e/OlJXv/jdfZn7+fZQc/i7eJt79CalM7dncB778HnxglkvbyYnHffJX/tWgLuuQffm25EGFrmeB9K/Qgh8A32wDfYgx5DI6gsN1vbZmiNQP/ep0134OnnSnC0F8Ex3gRHe+MfbkSvVyUaSuPk5eXpH3rooYwRI0YUlZWVifj4+Ojp06eHr1ix4nBt+990003RvXv3Lvrmm2+SVq9e7X3rrbfGJCQk7G/IfCG2UllZieEifI/Wt/GmBW2G07PYYuRNIYQerWjtmJRylBAiClgJ+AG/AbdKKc8541FLqQqpSUrJJ4mf8PTPT9PGvQ0vXfkSnfw72Tusi6bs4EFOPPscJT/9hHNUFEGzZ2MccqVq4KmcRUpJ3okSDh/IISMlj+Mp+RTna18ZTgYdQZFeBEd7awlHtDdunmoU2CotqSpk7NixURs2bPBzdnaWOp1Ozpo1K/3JJ5880RTnWr58uc+CBQtCExMT/zzztT/++MOld+/enU+cOLHX19fXAnDppZd2uOmmm3LmzJmTdeb+8fHxke7u7pYjR44479692zMmJqb0o48+OtS5c+dygK+++spj5syZEWlpaS7t2rUrf+mllw4PGzasGLQSmqVLl/49duzYQoBZs2aFpqSkuKxfv/5QQkKCc8eOHbu++OKLac8991xIWFhYxQ8//JA4ceLEyO+//97bbDbTrl278s2bNye1bdv2QhOeOn9n6lti8QN1JBY2MgM4CHhZnz8LvCSlXCmEeA24A3i1Cc/vkIQQ3NjhRjr6dWTmdzO5dfOt/Pey/zI6ZrS9Q7soXDt1IuLddyja+h2Zzz/P0X/9C/d+/Wgzdw6ul1xi7/AUB1KzNKP71W2RUlKUW87x1HxtScln71eHsVi0rzHvQDct0YjRkg2/UCM6nUpYG+Kb9w62zTlW1KTd2PzCjCVX39bpyPn2W7du3aGwsDBjzRvtmZKSkpx79epV5xfI888/f/juu+8+5wyhAN9//71n+/bty2p7be/evW7h4eHlVUkFQOfOnUsPHDhQ54hwGzZs8Fu3bl3iwIEDS+Lj46Pmzp0btmnTptQTJ07o4+Pj455++unDd955Z84777zjFx8fH5eYmLgvODi4XvXkP/zwgzEhIeGATqeTS5cu9S8sLNQfOXLkDzc3N8vOnTvdPTw8bDoldb0SCynllbY8aU1CiHDgOmABMMvaWPQqYJJ1l+XAfFphYlGlW2A3Vo1axZwf5vDwtof5I+uPFj0UeE1CCDyvGoJx0EByV63i5P+WcCh+PN5jxxJ4/wwMbWw6K7LSQggh8PRzxdPPlbje2u+IqcJMZlphdbJx+M9sEnYdB8DgoqdNlBdB7bzwC/XAP8yIbxt31Si0BYqLi6soLCzc25hjfPrpp16rV6/237Zt28HaXi8sLNR5enqedtP39vY2p6en1/mlfc011+QOGTKkBOCWW27JmTdvXjjA6tWrvdu1a1c+ffr0HIC77ror59VXXw1atWqVz3333Vev6c4XLFiQ7uXlZQEwGAwyNzfX6c8//3Tp169f6aBBg0rqd9X1V6/EQghhBhZJKR+osW0OcIuUsrHTdC4C5gBVDVv8gTwpZVWxzFEgrJHnaPb83fx5fdjrLP5tMe8eeJeDOQd5cfCLtPFoHTdWYTDgd8steF9/PSdff53c996n4Isv8P/HP/C/Yyo6Dw97h6g4OCdnPaFxPoTG+QDW2XhPlp0q1Ug9vVRD6AQ+QW74hRrxD/PAP9SIX6gHXoFuqnSjhvqUJLQk33zzjcfUqVOj3n///ZRu3bqV17aPp6enpaio6LRmAgUFBTqj0VhnCUObNm0qqx57eHhYSkpK9ADp6enO4eHhp50nPDy84tixY/X+zzImJqb62NOmTcs5cuSI86RJk6ILCwv148aNy3n55ZePubi42KxWor5VIYIaU6VbxQFdGnNyIcQoIFNK+asQ4soa5zpTrRcshLgTuBMgIqLlj/ngpHNiVu9ZdAnown+2/4cJGyfw1KCnGBg20N6hXTR6Ly/azJ6N78SJZC5cyMlXXiHvk08IvH8G3mPHIvQOMdmu0gwIIfAOdMM70I0O/YIBMJss5GWWkHOsmOz0InLSi8k6UkjKnszqbyEngw7fEA/8Qz2qkw6/UCMePs6q/U8zkJSU5Ny9e/fOdb2+cOHCtGnTptVaFbJ9+3a3CRMmxC5duvTvMWPG1FrVAtCjR4/So0ePuuTm5uqqqkMOHDjgfuONN563iuVMoaGhFRs2bDitV+axY8echw8fng/g5uZmKS4uri5aO378+Fn3dZ1OV30PdXFxkQsXLsxYuHBhRkJCgvPIkSPjXnnllbKZM2eevNDY6nLOxEIIkVrj6VQhxNiqOIFwIK+R5x8AjBZCjARc0dpYLAJ8hBBO1lKLcLTJz84ipXwDeAO0xpuNjKXZGB45nFjfWB78/kGmfT2NO7rcwfSe0zHoWn7VSBXn8HDCX3qJkttuI/OZZ8l4+BFy3v+ANnNm43H55fYOT2mm9E46/EON+IcaieNUaWBluZmcjGJy0ovITi8m51gRhw/m8NdPx6v3MbjqMfq6YvR1sS41Hvu4YvRzwdm1vv/LKQ0VEBBQmZyc7ALUeuOPi4urKCkp2XOhx929e7fr6NGj2z/77LOHJ02alH+ufbt161besWPHkrlz54YuWrTo2OrVq70TEhLcJk+enHuh542Pj8+fN29exGuvveZ3xx135Cxbtsw3OTnZdcKECfkAl1xyScnKlSv9xo8fX7Br1y63zZs3+15xxRV19tPfuHGjZ5s2bUy9evUq9fHxMTs5OUm9Xm/T++f5fssjrWuJdtP3OuP1txpzcuvMqA8BWEssHpRS3iKE+AQYj9YzZAraHCVKDdHe0awYuYJndz/L2/vf5tcTv/L84OcJ9gi2d2gXlXvPnrRb+RGFmzeTufBFDk+9A+PgwQTNma1mUFVsxuCip02kF20iT/8KLCuqJCejiOxjxeSdKKEor5yinDKyjxZRUlhxVlmrs5tTjWTDBaOflnx4+Ljg4mbA2U2Ps6sTBlc9Bhe9KgFpgNmzZx+fPXt228ceeyx85syZGY8//rhNeoU8++yzwbm5uU4zZsyInDFjRiRopQnJyckHACZNmhQBUNX9dNWqVam33nprlJ+fX8/g4OCK999/P6UhXU2Dg4PNq1evTp41a1bb2bNnR0RERJSvXr06OSQkxATwzDPPHLv55pujfX19e/Tt27dwzJgxObm5uXXe29PT0w333ntvuxMnThjc3d0t119/fc60adPq1Vajvs7Z3VQIMRitauJbYDWw1PqSGTgipUyzWSCnEotRQohoTnU33QNMllLWWpdVpSV2N62vzYc2M3/HfAx6A08OeJIr215p75DswlJeTu4HH3Dy1dfUDKqK3ZlNForzyrVkI7eMotxybckpozivnMLcckoL6u5FL4SWiBhctWTD2VWP3qBDr9ehc9IhhNYORAjBJQNDiLjEv0FxtqTupspF1ehJyAYDR6WUKbaMypZac2IBkFaQxuzvZ3Mw5yC3XXIb9/e6v1X0GqmNmkFVaS6qko/ivHLKS01UlpmpKDNRUWpdl5moKDVRUWamotSE2WTBbJJYzBakRWuAKi2SPtdFEdenYQ25VWKhNFCjEwsn4FagJ2Cs8ZKUUt7R6PBsoLUnFgDl5nJe2P0CKxNW0jWgK89d8Rzhnq13vo3TZlANDSVw5ky8Rl6rGngqSg0qsVAaqNGJxfucGleiZqWftMXIm7agEotTvkr7ike3PwrA4wMeZ2i7oXaOyL5qzqDqEhdLwPR78Bw+DKFTYxQoikoslAZqdGJRgNbQcxVaT5DqN0kpZ9ogwEZTicXpjhQeYc73c9ifvZ+JHSfyYO8Hcda33mGMpcVC4RdfkLVkKRWpqbh07EjgvfdgvOoq1UBOadVUYqE0UKMTiwRgu5Ryqi2jsiWVWJyt0lzJi7++yAcHP6CTXydeGPwCEV4tf7yPc5FmMwWffUbW0qVUph3GtUsXAu+7F49Bg1SCobRKKrFQGqjh06ZbLQauF0KME0JECyEiqhbbxKc0BYPewNy+c1k8ZDHHio5x46Yb+eLQF/YOy66EXo/36NHEfPYZIQuexJyby5E77yJt4iSKd+ygPom2oiiKUrfGzm4qpZQOMeKLKrE4t/SidGb/MJs/sv5gQvsJzOkzB1enOufDaTVkRQV5az/l5GuvYTp+HPfevQmccR/uffrYOzRFuShUiYXSQI0usYBTw3rXXFTrt2Yi1BjKsmuW8Y8u/+CTxE+Y9PkkEnIS7B2W3QlnZ3xvvomYLV/S5pFHqEhLI+3W2zg8dSoley54cD5FUZRWr16JgZRSV9fS1AEqtmPQGZh16SxeufoVckpzmPjZRN7d/y5mS71m3m3RdM7O+E2+hZivthA0by5lCYmkTZzE4TvvpHTffnuHpyhKC9a3b98OL774YoC947CVeicGQgiDEGK4EGKaEMLF2sbCvSmDU5rGoPBBrB2zlivCr+DFX19k6pdTOVZ0zN5hOQSdqyv+t99O7FdbCHxgFmW//8HfEyZw5F/TKTtY6wzJiqI4kLCwsK7r1q3zPP+eDTd+/PhIIcSl+/fvr3PUvR07drh17ty5k5ubW8/OnTt32rFjh1tTxuRI6pVYCCHaAnuBzcD/AB8gGXii6UJTmpKfqx8vXfkSTw54koTcBOI3xPNp0qeq8aKVzt2dgP/7P2K++ZrAGfdR8ssvHLphHEfvm0F5UpK9w1MUxU6+/PJL499//33OYXzLysrE+PHjY2+88cbsnJycvRMnTsweP358bFlZmV27nlVWVp5/Jxuob4nFS0An4CRag88TwA/ANU0VmNL0hBCMiR3DmtFr6OjXkf/u+C8zts4gu9Sm89E0a3qjkYBp04j9+isC/vUvirdvJ3X0GI498CDlqYfsHZ6iKDWMHTs2KiMjw/nmm2+Oc3d37/nII480bJzzOlRWVnL//fdHLF269PC59vv88889TSaT+M9//pPp5uYmH3nkkUwpJZs2baqzJCUtLc25V69eHT08PHoOGDAgLiMjo7pjxIcffugdGxvb2dPTs0ffvn07/Pbbb9Ut788sOYmPj4+87777QgE2bdrk2aZNm24PP/xwcEBAQPcJEyZEZWRkOA0ZMiTW09Ozh7e3d49LL720g9ls2+rw+vbouALYBKQA91m3JQP9bBqNPZTlQ+IW6DbB3pHYTZgxjHdGvMP7f77Py7+9zLgN43js8sda7WRmtdF7eRF437343jqZnHeXkfPBBxRs3oz39dcTMP1fOEeontdK6/Tlq4vanjyS1qTV4gFt25WMmHb/kfPtt27dukNhYWHGpUuX/j127Nhap01PSkpy7tWr1yV1HeP5558/fPfdd+fU9trjjz/epn///oX9+vUrPVcc+/btc+3YsWOprsbovh07dizdt2+f2/jx42ud0nzt2rV+n332WVJ0dHTFVVdd1f6JJ55o88orrxz7448/XP75z39Gr1ixImXkyJGFTzzxRNDYsWNjExMTD7i6up63iDk7O9uQk5OjP3LkyB9ms5m5c+eGhoSEVJw8efJ3gK1bt3rYegyf+pZYSODMafgigSKbRmMPO/4Ha/8JCa17fAed0DGl8xRWjlpJoFsg9357L/N3zKe4stjeoTkUJ19fgmbNJParLfhNmULBF1+Qcu1IMv7zHyqPqXYqiuLo4uLiKgoLC/fWtdSVVCQnJxuWL18e+MILL6Sf7xxFRUV6Ly+v04oBPD09zYWFhXVOgTFx4sTsbt26lRuNRjlu3Licffv2uQO8//77fkOGDMm/4YYbClxcXORjjz12oqysTPf1118b6zpWTUIIuXDhwnQ3NzdpNBqlwWCQJ06cMCQlJTm7uLjIa665pkhn4+kN6lti8QswCq2UAiHEJ8AwtFKM5m3Qg5D4Bay7G+76EXza2jsiu2rv254V161g6d6lvLv/XXZl7OKpQU/RM6invUNzKE7+/rSZOwe/f9xO9ptvkbdyJXnr1uMzPp6Au+/G0MamJbCK4rDqU5LQEkyfPj1izpw5Gf7+/uetNzAajebCwsLT7tZFRUV6T0/POt8bHBxc3QDC3d3dUlJSogNIT083tG3btvofe71eT0hISMWRI0fqNX21r6+vyd3dvbpk49FHHz0+Z86c0GuuuaY9wG233Zb11FNPHa/PseqrvmnKHKAEuARt/Ip4IB942JbB2IXBFSYsB7MJVv8DzBencYsjc9Y7M/PSmSy7ZhkSye1f3M6iXxdRqX42ZzEEBRH88L+J+WoLPvHjyFu9hpRhwzn+1FOYsrLsHZ6iKGdISkpydnd371nX8uqrr/rV9r4dO3Z4Pvroo+EBAQHdAwICugMMGjSo42uvvXbW/l27di3766+/3C0WS/W2v/76y61r167nrEKpTWhoaOWRI0eqJ3qyWCxkZGQ4t23bthLA1dXVUlxcXH0vz8zMPC3hOLOaw9fX1/Lmm28ePXr06L5169Ylvfbaa23Wr19v01409R3H4gBa4805wCvWdWcpZcvo4O8fA2P+B0d3w9fz7R2Nw+jVphdrRq9hbOxY3t7/NpM+n0RSruoRURtDcDAh8+cTs3kzXqOvJ/fDFSQPG86J557HlFNryaqiKE0gICCgMjk5uc5eG3FxcRUlJSV76lqmTZtW6x/sX3/9tX/Pnj0HfvvttwO//fbbAYDVq1cnT548OffMfUeOHFmo1+vlggULgkpLS8VTTz0VCDBq1Kha232cy+TJk3O2bt3qvX79es/y8nIxf/78Ns7OznLo0KFFAJ06dSpdvny5n8lkYvXq1V67d+8+Z5Lw0Ucfee/fv9/FYrHg4+Nj1uv1Uq+37STl9e1u6gqUAgullPdIKV8Aiq3bW4bON0Cf/4OdS+DPDfaOxmF4GDx47PLHWDxkMZklmdy86WaWH1iuBtWqg3N4GKFPPknM55/hNWIEOcuWkTx0GMefXEDF4XM2JFcUxQZmz559fOHChSGenp49/vvf/9qsTjIsLMwUERFRvQC0adPGZDQaJcAVV1wRN2/evGAAV1dX+cknnySvXLnS39fXt+cHH3wQ8MknnyTXp7Hlmbp3717++uuvH5o1a1ZEQEBA982bN/usW7cuqepYixYtOrxlyxYfb2/vnh988IH/sGHDzkp0akpMTHQZMWJEe2vvk0633357VkMSnnOp71whOwF/oIOUUgqtbOUAkCulHGDLgBrKJnOFmMrh3WshKxHu/A4CYm0RWouRXZrNYzsfY+uRrXQN6Mqjlz1KB78O9g7LoZWnppL9+hvkf/45mEx4Dr0avylTcLv0UjWbquIQ1FwhSgM1eq6QLsD30pqFWNc/At0aH5sDcXLR2lvoDbDqVqhQPSJq8nfz5+UhL/PsoGc5VnSMmzbdxEu/vkSp6YKrDVsNl+hoQp99hthvvsb/rjsp2f0LaZNv5e8JN5K/6TPkRRqwRlEU5WKpb2JRCpz573t76/aWxactjH8bMg/CxhmgRqI8jRCCkdEj2TB2A2Nix/DO/ncYt34cO9J32Ds0h2YICiLo/vuJ3fotwfPnYykuJv3BB0kexiqClAAAIABJREFUNpzst97CnJ9v7xAVRVFsor5VIVuAq4H1wC60gbHGAF9LKUc0aYT1ZPNp079/HrY+CcMXwOX32O64Lczu47t5fOfj/F3wN9dHX8+DfR7Ez7XWRtVKDdJiofjHH8letoySnT8h3N3xueEG/G67Fed27ewdntKKqKoQpYHq/J2pb2LRH9gKuKANliWAcmCwlPJnGwXZKDZPLCwW+GQK/LUJbvkEYofa7tgtTLm5nDf/eJO397+N0WDkwd4PMjpmtGpDUE9lf/1FzvL3KNi0CWkyYbzqKvxvn4Jb797qZ6g0OZVYKA3UuMQCQAjRGZiGNuLmIeA1azdUh2DzxAK0NhZvj4C8w/B/36rGnOeRkpfCYzsfY0/mHvoF9+O/l/2XCC811HV9mbKyyP3oI3JXfIQ5Lw/XSy7B7x+34zViBMLZ+fwHUJQGUImF0kCNTywcXZMkFqAlFW9cCW5+8M+vwc3H9udoQSzSwurE1bz060tUWiq5u/vdTOk8BYOuXoPEKYClrIz89RvIWb6citRUnIKC8J08Gd8bJ6D3Ub9/im2pxEJpoEZXhXgDs4CeQM3xyaWU8upGh2cDTZZYAPy9Hd4bDdFXwqRVoLPtYCItUWZJJs/8/AxfpX1FrE8s8y+fT/fA7vYOq1mRFgvF27aRs2w5xTt2INzc8LlhLH633YZzZKS9w1NaCJVYKA3U6MTiM7Qp0s+s8JVSSoe4yzZpYgHwy7uw6X64/F4Y/mTTnaeF2Xp4Kwt2LSCzJJObOtzEjF4zMDrXa+4cpYayhERyli+nYONGrR3GkCH4TZmCe98+qh2G0igqsVAaqNGJRRFa19KlQB5aA04ApJQv2yDARmvyxALgswdh95sw9jXoMbFpz9WCFFcWs2TPEj48+CGB7oHM6DWD66KuQ69Kfi6Y6eRJcld8RO5HH2HOzcXlkk74T5mC17XXqnYYSoOoxMKxJSQkOHfs2LFrRUXFrwaDQ1UpN3qArERgk5RyvpRykZTy5arFNvE1E9c8DZGDtPEtDv9k72iaDQ+DB3P7zuXDkR8S4BbAw9seZvzG8Xx7+FtaShufi8UpIIDA++7VxsN44nFkeQXpc+eRfPVQTr7+Bua8PHuHqCh2FRYW1nXdunU2nVQLIC0tzXDVVVfFBgUFdRNCXJqQkHBaJl9aWiomTJgQaTQaewYEBHSfP3/+OYcTf+yxx4ICAgK6e3p69pgwYUJkaWlpiyl6rG9i8TEwQQgxSwhxlRDiiqqlKYNzOHoD3PieNojWipsgK8HeETUrXQO7svK6lSwcvBCTxcSMrTOYvHkyu4/vtndozY7O1RXfCROI3rSRtm++gUv79mS99BJJVw4h47HHKD90yN4hKkqLotPp5PDhw/NXrFiRUtvrDz74YGhqaqrLoUOH/tiyZUvCkiVLglevXu1V275r1qzxWrx4ccgXX3yRkJqaui8tLc3lgQceCG3aKzg/k8lkk+PUN7F4GnADnge+QhvTYivwrU2iaE7c/WDyGtA7wwfxUJBh74iaFSEEwyOH8+mYT5l/2XyOFx9n6pdTufuru/kz+097h9fsCCEwDhpExNtvEbVhPV6jriN/zVpSrx3JkbunUfTjj0gbfVkoiqMbO3ZsVEZGhvPNN98c5+7u3vORRx6x2SRkbdu2Nc2bNy9r8ODBtc718Mknn/g//PDDGYGBgeZevXqV3XLLLVnLli0LqG3fZcuW+U+cOPFk7969ywIDA80PP/xw+qpVq2rdt8prr73mHxIS0tXX17f73Llzg6u2l5aWiqlTp7YNCgrqFhQU1G3q1Kltq0o/Fi9e7H/ppZeeNqGTEOLS/fv3uwDEx8dH3nLLLRGDBw+OdXNz67lp0ybPjz/+2DsmJqazh4dHz6CgoG4NmcjNqZ77HaZGu4pWzzdSGzRr2XXw4QT4x+fgWmtiqtTBSedEfPt4rou+jo8TPubNfW9y06abGBE5gnt63EOkd6S9Q2x2XNu3J/TJJwmaOZPcj1aSu2IFRd99hz4wAO/rRuE9dgyuHTvaO0ylhclZndi28nixe1OewxDsUeI3vv2R8+23bt26Q2FhYcalS5f+PXbs2Fpn7ExKSnLu1avXJXUd4/nnnz9899131zp1el2ysrL0WVlZhj59+pRUbevRo0fp5s2ba+0fnpiY6DZ69Ojqesu+ffuWZmdnOx0/flwfHBxc69TR27dvNyYlJe3ft2+f6+DBgzvddNNNeb169Sp76KGHQn799VePPXv2/CmE4LrrroudN29eyMsvv5xen9g3bNjgt2bNmqSrr746uby8XLRr167r+++/n3rNNdcUZWVl6RMSEuqcgr4u9SqxkFJGSimjalsu9IQ1CSHaCiG2CiEOCiEOCCFmWLf7CSG+EkIkWde+jTlPkwjtoVWLZB2EjyeDqcLeETVLrk6uTOk8hc3jNnNXt7v44egPjF0/lvk7tNIM5cI5+fsTeM90Yr/bStj/FuPWvTs5H37IobE3kDp6DNlvv03liUx7h6kodhEXF1dRWFi4t67lQpMKgPz8fB2Av79/dVLg4+NjLi4urrWFeklJic7X17d6Xz8/P7P1OHW2aF+wYEG60WiUl112WWmHDh1Kf/nlFzeANWvW+P373//OCAsLM4WGhpoeeeSR9NWrV/vXN/ahQ4fmDR8+vFiv1+Pu7i6dnJzkvn37XHNycnSBgYHmgQMHlpz/KKerb4lF1VgWE4B2wN/AGillY1uKmYAHpJS/CSE8gV+FEF8BtwPfSCmfEULMA+YBcxt5LtuLvRpGL4F1d8P66XDD66Crb+2SUpOnsyf39LyHmzvezFv73uLjhI/ZmLKRSZ0mcUeXO/BxVQNDXSidszNew4bhNWwYptxcCjZvpmD9BjKff4HMhS/i0b8/3mPH4Dl0KDr3Jv2HU2nB6lOS0NJ5e3tbAHJzc/Xu7u4m0JINDw+PWksf3N3dLXl5edVJRG5urs56nFr3B4iIiKieCtnNzc1SVFSkB8jKynKOiYkpr3otOjq6IjMzs97dR8LDw0+bYnnlypUpjz/+eMgTTzwR3qFDh9Knn3766NChQy9oqu963QWFEO2BP4HXgX8DbwAHhBAdzvnG85BSZkgpf7M+LgQOAmFoE5wtt+62HBjbmPM0qR4T4er/wr5V8M18e0fT7AW4BTCv7zw23bCJa6KuYfmB5Vy79lpe//11SiovOHFWrJx8ffGbNInIj1cSvflzAu6+i4q0NNLnzCVx4CDS586jeMcOpLnO7zVFaRGSkpKc3d3de9a1vPrqqxc8i2JgYKA5MDCw8ueff67O0Pfu3evevn37str2b9++fenvv/9eve/PP//s7u/vb6qrGuQ8565ISUmprq44dOiQc1BQUCWA0Wi0lJaWVt/nDx8+fFZhghDitGYOgwcPLvnmm29SsrKyfh81alTu5MmTYy40pvr+e/0CEAL8jtZDZK/1+XMXesK6CCEi0Ub23AW0kVJmgJZ8AEF1vOdOIcQvQohfsrKybBXKhRs4C/r8E7a/DNsW2S+OFiTMGMaCgQtYM3oNfYL7sGTvEq5dey0fHvyQCrOqdmoMl6goAu+7j5ivttDug/fxvu46Cr/9lsNT7yD5qqvJfOEFyhIT7R2mojRIQEBAZXJycp3tAuLi4ipKSkr21LVMmzatzqqQkpISUXWjLisrEyUlJdVdRMePH5/91FNPhWRlZen37Nnj+uGHHwbcfvvtJ2s7zpQpU7I/+uijgF9//dU1KytL/9RTT4XceOONte57PjfccEPOM888E5Kenu6UkZHhtGDBgpD4+PhsgN69e5ckJye77dixw62kpETMmzfvnD1PysrKxKuvvuqXnZ2td3FxkV5eXha9Xn/B7Svrm1j0B7ZIKXtJKSdJKS8FvgQuu9AT1kYIYQTWAPdLKQvq+z4p5RtSyt5Syt6BgYG2CKVhhIBrn4Mu8fD1o7Bzqf1iaWHifONYfNViPhj5ATE+MTzz8zOMXjeajSkbMVvUf9eNIXQ63Hv3JuSJx4nb9iNhi17CtVMnspct59DoMaSOG0f2smWY7Jm0K8oFmj179vGFCxeGeHp69mhIj4Zz8fDw6OXt7d0ToEePHl08PDx6Vb22cOHC9MjIyPKoqKhuQ4cO7TB9+vQT48ePL4BTpSRJSUnOAOPHjy+45557jg8fPrxDVFRUt/Dw8IqFCxfWq7HlmZ555pmM7t27F3fv3v2Sbt26XdK1a9eSZ555JgOgW7du5TNnzky/7rrr2kdHR3cdOHBg0fmOt2LFCv+oqKiuRqOx59tvvx34zjvvXHDf9fqOvJkF7JRSjq6xbT1wuZSyUXd0IYQB2AR8KaV80botAbhSSpkhhAgBvpNSnrPa5aKMvHk+ZhOsmQp/rodrn4d+d9o3nhZGSsnO9J0s+m0RB3MOEusTy70972VI2yFqWGsbMmVnU/DZ5+Rv2EDZ/v2g1+Mx4HK8R4/B8+qr0Lm52TtExYbUyJtKAzV6SO8vgGHAd2jtIDoCQ9BKMa5taFRCuxssB3KklPfX2P48kF2j8aaflHLOuY7lEIkFgLkSPrkd/toE170Ife6wd0QtjkVa2JK2hSV7lpBWkEa3wG7c3+t++gT3sXdoLU55Sgr56zeQv3EjpowMdB4eeI4YgfeYMbj36Y1QjZWbPZVYKA3U6MSiG/A94I02noUA8tFKFX5vaFRCiIHAj8A+wGLd/G+0dhargAi0MTQmSCnP2QXIYRIL0LqerroVEr+A6xfDpVPsHVGLVGmpZH3yel79/VUySzIZEDqA+3rdxyX+dXZRVxpIWiyU/Lyb/A0bKPziCywlJTiFhuA96nq8x4zGJeaC23cpDkIlFkoDNS6xABBCBAO3ApFo3U0/qGpg6QgcKrEAMJXDykmQ/A2MfQV6TLJ3RC1WmamMlX+t5K39b5Ffnq8G2WpiltJSCr/5lvz16ynevh0sFly7dMF7zBi8rhuJk98FN6pX7EglFkoDNSyxEEI4Ae5AiZTSdL7t9uRwiQVAZSl8dDOkfg+j/we9brV3RC1aYUUhyw4s4/0/36fCXMHY2LHc3f1ugj2Cz/9mpUEqMzOr22OUHzwITk4YBw3C+/pReAwciN5LjUjr6FRioTRQgxOLxcBdQEcp5aEa2yOABOANKeUMGwbaYA6ZWABUlMDHt0DKtzB0Pgy4X+tFojSZk6UnefOPN1mVuAq90DOp4ySmdpmqBtlqYmUJieRvWE/Bxk2YMjNBr8etWzc8Bg7AOGAArl27IvR1Diyo2EkDE4tv0KrGldYpH7i6rhfPl1gkAylSyhG1vLYJ6CCljLNFlI3lsIkFaG0u1k2D/auh/3QY/qQaofMiOFp4lFd/f5WNKRvxMHgwqdMk4uPiCTXafRLBFk2azZTu3UvRtm0Ub9uu9SyREp23Nx6XXYZx4AA8Bg7EEKxKkhxBAxOLy4EnUMlFa5QP/AfYUdcO50ssSoH3pJR31fLa68BtUkqH6Hvm0IkFgMUCXz4Eu16DbjfBmKXaNOxKk0vKTWLJniVsPbIVgMtCL2Nc3DiGtB2Cs97ZztG1fKbcXIp37KB423aKt22rHhfDOTYG44CBeAwciHuf3uhcXe0caevUwMRCUep0vsTiJHBYStnrjO0C+A1oK6U851SvF4vDJxYAUsKPC+HbJyB2GNy4HJw97B1Vq5FelM665HV8mvwpx4uP4+viy/Ux1zMubhwxPqpXw8UgpaQ8MYnibdso3r6Nkl9+RVZUIFxccO/dG4+BAzEOHIBzbKwam+QiUYmFYmvnSyw2AiOB94EFQBraJGSPAJOBz6WU11+EOM+rWSQWVX5dBptmQtilMPFj8Kj3RHSKDZgtZnZm7GRt0lq2Ht6KSZroEdiDcXHjGBE5AneDmpDrYrGUllKye3d1tUlFaioATsHBeAy4HOPAgXhcdhl6H9U+pqmoxEKxtfMlFgPQBsWqrUGABRgipdzWNKFdmGaVWAAc3Air7wB3f4h/EyIH2juiVim7NJtNqZtYk7SGQ/mHcHdy59qoa4mPi6dLQBf1X/NFVpmeXp1kFO/ciaWwEHQ6XLt2qa42cevWFeFU74mZlfNQiYVia+cdx0IIcTOwBKjZOT0XuFdKuaIJY7sgzS6xAEjfC6unQk4qXDEbBs8FvfrCtAcpJXuz9rImcQ1b0rZQaiolzjeO+Lh4RkWPwttFtVG72KTJROkf+yjeto2i7dso27cfLBZ0np54XHaZ1ttk4EAMoaoxbmOoxEKxtfqOvOkGDECbZTQT2CGldKg5rJtlYgFQXgSb58LeD6BtPxj3Jvi2s3dUrVphRSGbD23m06RP2Z+9H2edM1e3u5r4uHj6BPdBJ1SPHnsw5+VRvHNndYmG6cQJAJyjo6u7tLr37avmMrlAKrFQbK3eI286umabWFTZt1prd4GA6xdBl3H2jkgBEnISWJu0lo2pGymsKCTMGMa4uHGMiRlDGw+bTpyoXAApJRUpKdVJRsnu3cjycoSzM66dO+PapQtuXbvg2qUrzpHt1Jwm56ASC8XWVGLhSHL/hjX/hKO7oeetcO2zqteIgygzlfHN4W9Ym7SWn4//jE7oGBQ2iHFx4xgUPgiDTnUdtidLWRklv/xK8fbtlP7+O2V//oksKwNAZzTi2rlzdaLh2qULhrBQ1X7GSiUWiq2pxMLRmCvhu6fhxxfBPxbGvwMh3ewdlVLD4YLDfJr8KeuS13Gy9CQBbgGMiRnDDXE30M5LVWM5AmkyUZ6SStn+fZTu20fZ/gOUJSRAZSUAej8/XLt0xs2aaLh17YJTYKCdo7YPlVgotqYSC0eV+j18eheUZMOwx6Hf3WoocAdjspj48eiPrE1ay4/HfsQszfQJ7sO4uHEMjRiKq5Ma8MmRWCoqKE9IoGz/fkr37ads3z7KU1K0wevQuri6duyIc1QUztFRuERF4Rwdjd7Xt0WXbqjEQrE1lVg4suJsWD8dEjdD3AhtllQPhxiPTDlDZkkmG1I2sDZpLUcKj+Dp7Ml1UdcR3z6ejn4d7R2eUgdLSQllBw9qpRr79lOelERFWhqyvLx6H523Ny6RkThHR+McFYVLdJSWfLRti3Bu/iO3qsRCsTWVWDg6KeHnN2HLI+DmAze8DjFD7B2VUgeLtPDL8V9Yk7SGr9O+psJSQSe/TlwRfgX9Q/rTPbA7BjWUu0OTZjOVGRlUHDpERWoq5YcOUXHobypSU6uHIwdAr8c5PLw64XCOiMApKAinwECcggJx8vdvFuNtqMRCsTWVWDQXx/dpY16cTIIBM+CqR9RcIw4uvzyfz1I/47NDn7H/5H4s0oKbkxu92vTispDL6B/SnzjfONV9tRkxFxVpCcehQ5SnplYnHBVpaciKitN3FgK9n5+WaNS1BGlrnYuLfS4IlVgotqcSi+akokSbyOzXZRDaC8a9AQEOMbmsch4FFQX8cvwXfsr4iV0Zu0jN14au9nP1o29wX/qH9Kd/aH/CjGF2jlRpCGk2Y8rKOrVkZmHKzDx9W1YWpuxsMJvPer/O0xOd0YjO3V1bPDxqfSwMBoTBSSsJ0elBJxA6He69e+MS17DvApVYKLamEovm6MA62HgflBdC1wkwcBYEqXr85uRE8Ql2Hd/FT+laopFZmglAuDGc/qH96R/Sn77BffF19bVzpIotSbMZc27u2QnHyWwsJSVYiou1dfVSjKVYeyxL6h6TMHj+o/jefHODYlKJhWJrKrForgpPwM7/we53oLIEOl0PVzwIId3tHZlygaSUHMo/xM6MnezK2MXu47spqixCIOjo11ErzQjpT882PXFzUqNKtlZSSqisRFZWIk0mpMWitcGyWLRSjQaOOKoSC8XWVGLR3BVnw65XYdfrUF4AccO1eUfa9rV3ZEoDmSwmDmQf0Eozju9iT+YeTBYTBp2BHkE96B/Sn34h/ejs3xknneM3DlQcm0osFFtTiUVLUZav9R7ZuRRKcyDqCi3BiBykxr9o5koqS9iTuYddGbv4KeMnDuYcBMBoMNInuE91iUaUd1SLHm9BaRoqsVBsTSUWLU1FMfzyLuz4HxQdh/C+WoIRN0wlGC1ETlkOPx//WUs00n/iaNFRAILcgqrbZ/QL6UeQe5CdI1WaA5VYKLamEouWqrJMmzF12yLIPwLB3bQ2GB2vBzUhU4tytPBodWnGroxd5JbnAhDtHU23wG7E+sQS7R1NjE8MwR7BqnurchqVWCi2phKLls5cCX+sgh8XQk4KBHaEQQ9A53GgV/XzLY1FWkjMTWRXxi52Zuzkr+y/yC7Lrn7dzcmNGO8Yon20RKPqcZgxTCUcrZRKLBRbU4lFa2Exw4FPtQQj80/wDIX2I7TGntGD1SyqLVheWR6p+amk5KeQkqctqXmp1V1cAVz1rkR5R2nJhk8M0d7RxPrEEmYMQ6/T2zF6pampxEKxNZVYtDYWizb3yN4VkPodVBSB3gUiB55KNPyi7B2lchEUVBSQmpeqJRv5WrKRkp/C8eLj1fs465yJ8o4i2kdLNKpKONp6tlU9UloIlVgotqYSi9bMVA6Hd0LiFkj6ErKTte3+caeSjIjLwKn5T7Sk1F9RRZFWwpGXUr1OyUshvTi9eh+DzkA7r3Za+w2faGK8tZKOCM8INRdKM6MSC8XWVGKhnJKdAklfaUnG39vAXAHOnhBzpTa7atxw8Gxj7ygVOympLOFQ/iFS8lNIzkuuLu04VnQMifY94iScaOfVjkjvSALdAglwCzht8Xfzx9/VXyUfDkQlFoqtqcRCqV15ERz6HhK/1JKNQut/qyE9IPZqaNMZAjqAfywYXO0bq2JXpaZS/s7/W0s2rCUcaQVpnCw9SUFFQa3v8XHxqU40AtwCCHC1Jh/uAac993bxVmNzNDGVWCi2phIL5fykhBP7rUnGFji6G6RFe03owKcdBHaAgPbWdQcIbA+u3vaNW7G7CnMF2aXZnCw9qS1l2vq0bdal3Fx+1vuddE74u/qfVepxWkmIq7bN3eBuhyts/lRiodiaSiyUC1dZqrXHyEqAk4mn1tnJWvVJFWOwlmAEdDiVeAS0B3d/1W5DOY2UkuLK4tMSkOzSbLJKsk57frL0JDllOViqEtsa3J3cq0s5PAweeDp74mHwwGgwYnQ2auuaj8/Y5ubk1iq73KrEQrE1h23WLYS4BngZ0ANvSSmfaYrzHMsr5fcjeVzTORidThW51ovBDYK7aktNZhPkpVkTjQTIStTWv6+EisLT93VyBRcvcPEEV68aj73r2O4FLt6nb3f2UKOJthBCCO1G72wk0jvynPuaLWZyy3NrLfXILs2moKKAosoiskqyKKosoqiyiOLK4vPHgMDd4I6L3gU3Jzdc9a64OLngqnfFzckNF70Lrk6nP3Z1csVVf2rtpHPCoDNg0BlOPdYbat9e47GTzgmd0KHX6dGhQye0RVUDKc2RQyYWQgg9sBQYBhwFdgshNkgp/7T1uVb+fJj/fZtM+zZG7r0qjpFdQ9DrBGaLpKTCVBUPOqF98Qih3csEp/7ga/7tN/RrQAjtiFXHsuUXSlWplJRwZvlU1Vm062rkOfVO4B+jLYysPre0SCjMQJ5MRGQnQ2muNuV7eQGUFyDKC6GsAIoya2wvRJwV7RnXJXRaolGVfLh4g6sXsvq5lzUJ0R5LgytC6LXqm6q1ruq57oznVa/r6thfr/3Qztq/9uOd9bM962ct6nxNnPVf9HmO1cJvRnqdvroapAMd6vUei7RQXFlMcWUxhRWFFFcWU1BRQHFFMYUVRRRXFlFYWURpZQll5nLKzKWUmcooN5dTZi6joKKAMnMZZSbrYn1caals0msVCPRCjxCn1lXbaz6e02cON8Td0KSxKEp9OWRiAfQFkqWUqQBCiJXAGMDmicX9Q9sTG2Rk8TdJ3PvRHh7beACTRZJfWomj1BJpiQw1vkhO3TtqJgvVCYR1e2PPeepctZ/31Lku5Lyh1uU858eCB2V4UoqnKMFIKV6iBE9K8BSlGK1rz8oSPItLte3kYxQZ1e/xpARnYb6wC1dsxiIdK8HxsC62nEHFDJQLQakQlOsEJgSVAkxCUCkElVQ9hkpx6vVKIbTtnNrXAtoiwIxACu34EoFZaK9J6/PqPzPrj7j0152gEgvFQThqYhEGHKnx/CjQ78ydhBB3AncCRERENOhEep1gTI8wRnUL5fN9GXz7VyZGFyd83Q14uhoQAixSIiVYJEjkaTfPmm1UGnozr7ohV3XZO61kQcrTbthV55ecUdpgfVYzIaCWhKCq1IVajlkVC7UkCmfuU9cxzzxvXfHZggQKrEtd9JZynE1FOJuLcTKXg7Sgw4KQ1oU61tb9kBZ059mv6nHd+58vuTn1iyPO+iU64/kZr5/945R1PG6ARr29EW9u1FvlqdLE+pQk1vgdrnmMuuKQNR7IM3Y47W+yRgwC7Yu2+su2lnMKqpKG04sVT/teAHTVpaai+m/SIiUBcSPrukJFuegcNbGo7XvgrD9zKeUbwBugNd5szAn1OsH13UO5vvv5/5tWFEVRFKV2jtoE+ijQtsbzcCC9jn0VRVEURXEQjppY7AbihBBRQghn4GZgg51jUhRFURTlPByyKkRKaRJC3AN8idbd9B0p5QE7h6UoiqIoynk4ZGIBIKX8HPjc3nEoiqIoilJ/jloVoiiKoihKM6QSC0VRFEVRbEYlFoqiKIqi2IxKLBRFURRFsZkWM7upECILSGvg2wOAkzYMx57UtTielnIdoK7FUTXmWtpJKQNtGYzSurWYxKIxhBC/tJRpg9W1OJ6Wch2grsVRtaRrUZo/VRWiKIqiKIrNqMRCURRFURSbUYmF5g17B2BD6locT0u5DlDX4qha0rUozZxqY6EoiqIois2oEgtFURRFUWxGJRaKoiiKothMq08shBDXCCEShBDJQoh59o6nMYQQfwsh9gkh9gohfrF3PBdCCPGOECJTCLG/xjY/IcRXQoiqtjmlAAAJe0lEQVQk69rXnjHWRx3XMV8Iccz6uewVQoy0Z4z1JYRoK4TYKoQ4KIQ4IISYYd3erD6Xc1xHs/tchBCuQoifhRC/W6/lMev2KCHELutn8rEQwtnesSqtV6tuYyGE0AOJwDDgKLAbmCil/NOugTWQEOJvoLeUstkN+iOEuAIoAt6TUnaxbnsOyJFSPmNN+nyllHPtGef51HEd84EiKeUL9oztQgkhQoAQKeVvQghP4FdgLHA7zehzOcd13Egz+1yEEALwkFIWCSEMwDZgBjALWCulXCmEeA34XUr5qj1jVVqv1l5i0RdIllKmSikrgJXAGDvH1CpJKX8Acs7YPAZYbn28HO1m4NDquI5mSUqZIaX8zfq4EDgIhNHMPpdzXEezIzVF1qcG6yKBq4DV1u0O/5koLVtrTyzCgCM1nh+lmX7hWElgixDiVyHEnfYOxgbaSCkzQLs5AEF2jqcx7hFC/GGtKnHoqoPaCCEigZ7ALprx53LGdUAz/FyEEHohxF4gE/gKSAHypJQm6y7N/XtMaeZae2IhatnWnOuGBkgpewHXAtOtxfKK/b0KxAA9gAxgoX3DuTBCCCOwBrhfSllg73gaqpbraJafi5TSLKXsAYSjlbp2qm23ixuVopzS2hOLo0DbGs/DgXQ7xdJoUsp06zoT+BTtS6c5O2GtH6+qJ8+0czwNIqU8Yb0ZWIA3aUafi7Uefw3woZRyrXVzs/tcaruO5vy5AEgp84DvgP6AjxDCyfpSs/4eU5q/1p5Y7AbirC2qnYGbgQ12jqlBhBAe1oZpCCE8gOHA/nO/y+FtAKZYH08B1tsxlgaruglb3UAz+VysDQXfBg5KKV+s8VKz+lzquo7m+LkIIQKFED7Wx27AULQ2I1uB8dbdHP4zUVq2Vt0rBMDaxWwRoAfekVIusHNIDSKEiEYrpQBwAlY0p2sRQnwEXIk2/fMJ4FFgHbAKiAAOAxOklA7dMLKO67gSrbhdAn8Dd1W1UXBkQoiBwI/APsBi3fxvtPYJzeZzOcd1TKSZfS5CiG5ojTP1aP8YrpJSPm79+18J+AF7gMlSynL7Raq0Zq0+sVAURVEUxXZae1WIoiiKoig2pBILRVEURVFsRiUWiqIoiqLYjEosFEVRFEWxGZVYKIqiKIpiMyqxUByGEOI7IYQUQtx+nv2kdYlswlj+tp7jyqY6x3nOv8x6/vn2OP+5CCGGWGN7oomO/3/1+T1QFMUxqcRCqbcaN9uq5aQQ4kshRG8bnWI18DLwp/V8863nWXbGfi9bl6YcXvod6zmONuE5mqunABOwpImO/x5wEnjMOgOxoijNiNP5d1GUs2wCDgGD0Ub47COE6GgdSrzBpJT1ulFJKe9vzHnqeY7Hm/oczZEQogfaENJbpJQnmuIcUspyIcSnwP+hzXuzqSnOoyhK01AlFkpDvC2lvA+42vrcF7gMtJEBhRBfWEszsoQQG4UQHareKIS4XwiRIoQos77+XdXrNatCrFUAj1rfNsW6/TvrfqdVhViHOX5LCHFYCFEghPhJCHFNjXNWVSu8Zo2nxDqjZY+6LvDMqpAasT0thPjBeoztQoh25zhGuBBiuRAizXq9B4UQfeoTcy3HOq1qRAgRWfVzqLFP1c9ljhDikBAiz/p4kBAiwfp8cY39q0qEVgsh3hNCFAkhkoUQQ+uKAxhlXX9/Rnz9hRBbrJ9pkfV63GvGKYS4VwhxXAhxQghxqxAi3nr9WUKIeWecp+r4o1AUpVlRiYXSIEIIHVqJRZWTQpt74XtgBPAT2tDCo4DvhBC+QohY4CXAC1iGNuVzBFBzzoYqP3FqauuDaNUSq+uIYwNwB1rx+XrgUuAzIcSAM3a/C60I/xDQFfjfBV20Zg5wxHquy4Ena9tJCOEOfAvcBpQB7wO5QOgFxtwQDwA7AW/gGbSf20+AC3BvLYlDPBCKNldGDFo1UF26WdcHqzYIITqjTYY1zLr9Y7QhzZ3PeO/9aJ9pENqkX0uAHwB/4CkhRPsa+1Ydv+c5r1RRFIejqkKUhvj0jOcb0W5kDwI+wHdSylEAQog9aPMxTECbrwG0mRfXAn9KKY/WVo8upfxCCNEf6Af8fI7qj95oRfNFwCApZbEQ4iTaTWw6sL3Gvp9LKW8QQgxBu+k35Kb1mpRyuhDiH2g34LqOMRKIQ5uOu6eUsgSqZ9m8kJgb4gEp5QdCiMuBdsByKeUcoU0bPs4a89c19j+AlhREAqlAWyFEgJTyZC3H9rWuC2tsm4aWtGyQUo6xXqcebQ4Onxr7TUFLcEqt+z8hpXxFCNEJ6AV0BxKt+1a1n/FFUZRmRZVYKA2xCW3itv+g3UDHWKeejrS+frDGvn9Z1+2klAfRqjfCgC+BI0KIv4BOjYil6pxHpJTFZ57zjH33WNd51rVHA8535jGMdewXZV3vq0oqAKSUlVxYzHU5V6PGqp9/VYwJ1nVVMnDmde+V2qRBeTW21XVdVft41thWda0/VW2oMR35aXFJKU1A1TWfKy6vM86nKEozoRILpSHellLOlFL+f3t3DOJEEIVx/P+JWqhgI9ddb3HWFjbmrEQ4EUsLQbTTwsrGUgttLARRmwMVEa4TTkWwEK1ECIidhcpBwEawVDyfxZvhYjA5k0xQ4fs1Q3aHmd0lMG/3vWUvRcTj2PiS3YfS7u3rW+srPpa72MsRsYdcQK+U/eeHzLNe2lH/0zrnfEk//DLnQN/vpZ3my3t/Osb70u5Tft4aAElbGe+Yq7oY1wV3YcTc65v8HjTOdXlT2v5gsJ7r/rpB0hZJmuK46vjdEX3M7B/kwMJaugd8ATqSHkp6Qj52/0Tm+eeBnqQV4AJQixWH3ZWulfawpOuSjv+mz2syb78LeCHpDnCOXCRvNDinST0C3pH1I11JtyQ9B44w2THXBfakpKsj+s1afUOjv77mJvAVOFqKXG+T6ZXdU8xTx1+dYgwz+wscWFgzEdEDOsBT4ABZS7AKdCLiM5k3f1X2nSELBh8wpAASWCFTJjuBs2XswTl/AEvAMlkUeIxchJci4mWrcxtXSX8cIos2d5D1BXNAb8JjvgvcB7aRBbHXZnoCQ0REl0x5LEqaK9veAgfJuo0F4AQZYH6bZA5J28lrsoYDC7P/jjaeYpuZbU7SIvCMTGtdnMH4p8m3Rk5FxHLr8c1sthxYmJmZWTNOhZiZmVkzDizMzMysGQcWZmZm1owDCzMzM2vGgYWZmZk148DCzMzMmnFgYWZmZs38BMIjqQMnlw6qAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "%matplotlib inline\n", "\n", "# Define an array of x values to estimate the function\n", "positions = L * np.array([0, 0.01, 0.02, 0.04, 0.06, 0.08, 0.1,\n", " 0.15, 0.2, 0.25, 0.3, 0.4, 0.5, 0.6, \n", " 0.7, 0.8, 0.85, 0.9, 0.95, 0.98, 0.99, 1])\n", "\n", "# Define an array time points to estimate the function\n", "times = [0.00001, 0.1, 0.5, 1, 2, 4, 10]\n", "\n", "# Define a function to use to compute the value of the \"ith\" term\n", "# in the series of eigenfunctions that are summed in the solution\n", "def eigenfunction(Pe, B, x, t):\n", " return (B * (B * np.cos(B * x) + Pe/2 * np.sin(B * x)) /\n", " (B**2 + Pe**2/4 + Pe) / (B**2 + Pe**2/4) * np.exp(-1 * B**2 * t)\n", " )\n", "\n", "# Store the results at each x,t pair in a list of lists\n", "Cs = []\n", "\n", "# Estimate the concentration for each dimensionless time and x\n", "for t in times:\n", " tau = D * t / L**2\n", " Cs.append([])\n", " for p in positions:\n", " x = p / L\n", " \n", " # Get the eigenfunction values for all the eigenvalues\n", " series = eigenfunction(Pe, np.array(betas), x, tau)\n", " \n", " # Sum the series and convert the result to concentration at the point\n", " C = C0 * (1 - 2 * Pe * np.exp(Pe/2 * x - Pe**2/4 * tau) * series.sum())\n", " Cs[-1].append(C)\n", " \n", "# Plot the results\n", "fig, ax = plt.subplots()\n", "ax.set_xlabel('Position in column (cm)', size = 12, weight = 'bold')\n", "ax.set_ylabel('Concentration (mg/L)', size = 12, weight = 'bold')\n", "ax.set_title('Column Concentration Profiles', size = 14, weight = 'bold')\n", "for t,C in zip(times,Cs):\n", " ax.plot(positions,C, label = 't = {:.1f} hours'.format(t))\n", "leg = ax.legend(bbox_to_anchor = (1.02, 0.5), loc = 6, fontsize = 12)\n", "leg.get_frame().set_linewidth(0)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The figure shows the evolution in the concentration profile over time throughout the column. At the beginning, it is zero everywhere, as expected, while at the end, it has reached the influent concentration everywhere. Between these extremes, the concentration gradually breaks through the column.\n", "\n", "### Calculating the concentration in the effluent\n", "While it is interesting to plot the distribution of concentration throughout the column, the concentration at the outlet is what is of primary interest for tracer studies, since it can be compared to observed data. A high resolution time series of concentrations can easily be obtained by evaluating the function at the outlet ($x=L$) as shown below." ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYsAAAEaCAYAAADg2nttAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzs3Xd4VNXWwOHfEgIBErogPaAoIIYWQEQBQbBh73qVoiIqYkc+G8VykWu/CIoNFBQRC0W5KihSRUBROiK9SO89yfr+2CdhEibJpExOynqf5zwzc+o6k8msOXufvbeoKsYYY0x6TvE7AGOMMXmfJQtjjDEZsmRhjDEmQ5YsjDHGZMiShTHGmAxZsjDGGJMhSxaFmIisFREVkRF+x5IfiUg77/1TEWmX29v7SUSmeXFP8zsWkzssWeQzIlJcRB4WkTkisldEDovIXyLyvojU9zu+cBKREQFfrioiCSKyU0Sm5uUv24C41/odS15RmD/H+VVRvwMwoRORcsBUoIk36wDwF1AD6A4sApb5E12um4v7/J4DtAfOFZF6qrohrQ1EpJiqHsutAMMtv55PbnyORUSAIqoan539mBPsyiJ/GcKJf7CXgfKqGquq5XBfmIuSVhSR80XkO+9X21ERWSEiT4lIRFo7F5GYgF/tXQPmpyiuSrXeYyLypYgcEpHF3nGbiMivInJQRGaKyFkB+0r+lS0iN4jIcm+96YHrZURVz1XVOOAub1ZJoEWQ+PqIyNcicggY7C0/TUTeE5FNInJMRNaJyCARKR4Q5+3eOewQkeMistt7P1ukF5eI3O8dN1FE7vGuJrp4i2ulU+xUTUS+8t6LNSJyZ8A+A4ur7haRn0TkCHCft7yh9zfY4Z3PGhF5WUSiAvZxUrGRiPRP2m/AvGIi8l8R2SMiu0TkTRF5MfV6qc75Lu+Y+0Vkkoiclt57RIif4zTiO6noLnA9EblURJYCx4HHvXnxIlI5YB9J8/eKSAlvXnMv9l3e/8siEemWwXkULqpqUz6YgDK4fwAFFgKSzrrtAtbdDazwnivwacB6a715I7zXMQHrdQ1xvSPAatyvQwW2Atu9Yx7z5s0K2NcIb95xb/kyIDH1emmcV9K26r0uCgzw5iUCsUHiOwrsxX0BvQxUCDifA8Af3joKTAw41hDgsHceC73zVGAfcFrA+5x0nHZANy+OBOBOb52vvPcjKZZfvKlpqu0PAWu8WNXbR70gxznq7W8p0BuoD+wPOJ+l3rYKzAJO8fYxzZs3LeAc+we+n968wQHHWusd60CQ9ZL2d8h7n1YGbDc6hz7HweJL8Z6nXs97f1YD64FzcZ9/BXoF7GOBN2+49/q8gM/AVmBJwP4e9ft/P69MvgdgU4h/KGge8AH+bwbr/uyttx4o580bFLD9Od68tWQ/WXwHCO4XftK8d731nguYV8KbNyJg3hXevFdTr5fGeQVum3p6OmC9wPiWAWW9+UWAZ735O4Eq3vzWAeu39uadBZQM2OcZAeskJYLAL65huC/pBKBLGnGvTTU/cPtx3vsYGzCvZ5D1fgIiA85npDf/IFDLm98zyHs8jQySBe7q7LA3bzyu5CGKgEQQsG3S/hKAxt68L715/+TQ57h/kOMGvhftUq8HDApYtwjwjjd/pjevbsC653nzfvRe/wxEePOe4sSPg0i////zwmTFUPmHBDzXDNZt7j3+T1V3e88/CVgel2NRwTfeN83agHkTvcfVAfMqpdpur6omrbc0nfXSMheYj/tnBnhaRM4Lst5IVd0DoKoJQEtvfnlgs1fEMTNg/XO9xzLAeK9YIhFXpp6kapDj9MR9ufZR1ZEhnkOgUd77GPheVA6y3juqegSSzyfpbz1LVdd5z7P6tz4DiPSef6aqiap6AJiUzjaLVHWh9zwp9vT+hpn5HGfF68k7d+/PCO/leSJSA7jZe71SVWd7z5M+E22AY95n4nlvXjRwdhjizHesgjv/WAHE4/5m54uIJP0cTEdm/xkD1y8S8LxMOtskfVnHB5kXuL/ALwmAPQHP49NZLyhVPRdARCoC63C/ih8AZqda9Z9Ur5P2fwBX3JDaHq+s/zugLK746Xdc0UnSl0qRINsdwP0Kv0dEPlbVbaGcR+BxAVQ1XiT5LQj2XqQ+nyShfhZC/buG+tkJ9ndM72+Ymc9x8nwRKeJ9+acXM6r6T6rXc0RkBe5K8SZOJIsRQTbfDAS7QSIxvWMWFnZlkU+o6l5grPeyCfCiiCQnexFpIyLtvZfzvMdLxd15AnBrwO7mp3GYwC+40739tsN9aeZlSV9OaVbeB/jVe1TgX+oqys8FLsTVaXyB+2JJOufuqtoMeCiD/d4HbMIVc3wnIoFfaoe8x5ISkAlySNLfurWI1PKeB/tbJ/1tY0SkiIhEApem2tcqXHIEuE6cKKBzTgWbyc/xSZ9H4JosHDbpSu9RoAHuy/+jgOVJ7+FmoEPAZ+IK4HVV/T0LxyxwLFnkL71wv3IB+gK7RORPEdmBK2+N9Zb1w/16qwGs9n5ZPeEtG6OqiwhCVQ8Dc7yXj4nIT7gipTz3y0pEfhGRX3FfcCW82V+HsOkQ3K/HaGCp9/79hasI/RyXJFbj6gAA3heRP0PY9wbcl+8eoDEwKelOG2C593gqsNyLvUSQfWTFINxVTSlgiYgsAd7yls0GvvGeT/UeqwO/AYtxiS2Zqh4C/uu9vA73PqwBquVQrElC/Rz/xInP3o8iMhe4PQvH+8jbT9JdWj+o6qaA5U/jrhzjgC0i8ruIrMddxQ3KwvEKJEsW+YhX/3Ae7hdS0i/kM3F3w4wEvvfWm4b7pfw97m9cG1dJ+QxwRwaH6QrMwCWbarh/7DTbLvioJa68PgL3xXOfqo7KaCNV3YGrl3gP98u1PlAa9+vySWCr9z7fgCuDPwV319YVIex7EXA17s6a84EvxN2q/AHuimUv7u/VkuBFWZmmqsuAVri7ro56+98AvAJcrKpJX7YfAm8CO4CauErdN4Ls8mlcQt0LlMMl0A+8ZUeCrJ+VmEP9HC8HeuDqwyrgEvp9WTjeJuCHgFkjUi2fCVyAq5uJx119gEu0z2T2eAWVZFzsbYwpLLz2CEe84iK8K6B5uEreOaoa7CYCUwhYBbcxJlArYJSIzMP90m+OK76Jx912bAopK4YyxgRag2u0FgtchisumwC0UdUpfgZm/GXFUMYYYzJkVxbGGGMyVGDqLCpWrKgxMTF+h2GMMfnKggULdqjqqRmtV2CSRUxMDPPnp9XWzBhjTDAisi7jtawYyhhjTAgsWRhjjMmQJQtjjDEZsmRhjDEmQ5YsjDHGZChXkoWIfCAi20RkccC88iLyg4j85T2W8+aLN+7vKq8nyqa5EaMxxpi05daVxQjgklTz+gJTVbUurvvkvt78S3FdJ9fF9Tg5LJdiNMYYk4ZcSRaqOh3YlWr2VZwYlGQkrmvnpPkfqfMLUFZEquRGnMYYY4Lzs86isqpuAfAek8btrUbK8RM2ksbgKyLSQ0Tmi8j87du3hzVYc0L//u6xalUQOXkKXL55M0ycGHw9kRPLr/BGi7jiiuDrBS6fONFtl9Y+k5ZXrXoi3mDrBS63c7JzypFzGrqBzd8vpuqpx2D6dPrfsTr4OZU9CG++Sf9Os+nfYQYMHEjV6H3Bz6nRV3DrrVQtsYvNl97JxBYD0zzvcMq1jgRFJAaYpKoNvdd7VLVswPLdqlpORL4B/u0NSIKITAX6qOqC9PYfFxen1oI7dwT+0xqT56jCkSNw4ICb9u93jwcPnph38CAcOnTiMXA6fPjEY9J05MiJx6QpISFn4z7lFIiIgGLFoGhR9zwiIuXzIkVOzKtRA774ItuHFZEFqhqX0Xp+dvexVUSqqOoWr5gpabzdjbjhQJNUx42Na/IISxQmrOLjYfdu2LXLTbt3w549J6a9e93jvn3uedLj/v3u+b59mfsij4iAkiXdVKIElCrlHkuWhIoVITLSvQ58LF7cPS9e/MTrYFOxYu4xIuLE66TnSYkhad4pefvmVD+TxQSgC26M2y7A+ID5vURkDG74yb1JxVUmb0i6xDcmQ6rul/y2bbB1q3vctg22b3fTjh1u2rnTPe7a5b7s01O8OJQpk3I69VQoXRqio088RkWlfCxVyk1RUSeelyzpvqhNhnIlWYjIp0A7oKKIbAT64ZLEWBG5E1iPG/MY4FvcoCurgENAt9yI0YRui6VuAy4JbNzopk2b3C+IzZvdB2TLFvjnHzcdOhR8++hoqFDBfdGfeiqcdZZ7Xb78icdy5dxUtqx7LFPG/ZI3uS5XkoWq3pLGog5B1lXg/vBGZIxJl6r75b9mjZvWroV162D9ejdt2OCKflIrWxaqVHFTq1ZQuTKcdhpUquSeV6rkpqTiHZNvFJguyo0xmaTqrgRWrnTTqlXw11+werWbDh5MuX758lCzJtSpA23bugrW6tXdVK2aK58sWdKfczFhZ8nCmIIuPh7+/huWLIGlS2HZMli+HFasSJkQihd3ieCMM6BDB/e8dm2IiYFatVyxkSm0LFkYU5Ds3AkLF8Iff7hp0SKXII4ePbFOzZpQvz6cfz7Uqwdnngl167orhDx+R47xjyULk2n9+vkdgQHcnUPz5sH8+W767TdXn5CkShU45xzo1QsaNnRTvXrubiBjMsmShck0a2fhg4QEd5UwZw7Mng1z57r6hSR160Lr1i4xNGkCjRq5O4yMySGWLEymWTuLXHDsmLtamDYNpk93CWL/frescmV3p1G3btCyJTRr5m4pNSaMLFmYTLNeVcIgMdHVNUyZAlOnwsyZJ9onnH023Habu3Jo3dpVOIv4Gq4pfCxZmExbsOBE524mG7Zvh+++g8mT4Ycf3GuABg3gzjvhwgvhggtcmwRjfGbJwmTalVe6W/RNJqm6O5MmTHDT3LluXqVKcPHF0KkTXHSRq5g2Jo+xZGFMOKm6O5a++AK+/NI1fAOIi3O3lV1+OTRtaresmjzPkoUxOU3V1T+MGQOffea6yShaFNq3h0cfdQMjVAs6RIsxeZYlC2Nyyvr1MHo0jBrlipuKFoWOHd29xldd5TrCMyafsmRhTHYcOeKKlz780N3FpOpaRg8bBjfc4HpPNaYAsGRhTFasWAFvvw0ffeRaUteqBc8+C126uP6UjClgLFmYTNu0ye8IfJKQAJMmwZAhrj1ERARcfTX06OHqI6yS2hRglixMphW6dhb798MHH8Cbb7quu6tXh+efh7vucq2pjSkE7KeQybThw/2OIJds3QpPPeV6aX3oIdf+4fPP3WBATz1licIUKnZlYTJt4kS/IwizDRtg8GB4913XR9N118Fjj7l+mIwppOzKwmTaFVf4HUGYbNgAPXvC6ae7yut//ctVZH/+uSUKU+jZlYXJtEmT/I4gh23bBi+84BKEquuX6f/+zxU/GWMASxamMDtwAF59Ff7zHzh8GLp2hWeecbfBGmNSsGRhCp/ERBg5Ep58Ev75x9VJvPACnHWW35EZk2dZsjCFy+zZ8MADbgjSc891ra9btfI7KmPyvJAquEXkFBE5S0TOFZF6IlIk3IEZk6O2bnXFTK1buzqKTz5xicMShTEhSTdZiEgHERkL7AGWArOAJcBuEflcRNrnQowmj+nc2e8IMiEx0TUMqVfPJYi+fWHZMrjlFhttzphMSLMYSkR+BNoCAhwGFgP7gNLAGcB1wLUi8rOqWtIoRPJNO4sVK1wr65kzoW1bd7dTvXp+R2VMvpTelUUc8BbQGiitqo1U9QJVbYRLGK2Bod56phDJ8+0sEhLcHU6NGsGSJfD++/DTT5YojMmG9Cq4a6rqnmALVDUBmAPMEZFnwhKZybN69PA7gnT89RfccQf88ovr5G/oUBum1JgckOaVRVqJIomI3CAivTNazxQ8zZr5HUEQqvDOO9C4MSxf7uonvvzSEoUxOSQ73X08AryWU4GY/CPPjQi6axdce63rqqNVK1i0yCqwjclh1jeUyd9mzXJXE998Ay+/DN9/77oQN8bkKN+ThYg8LCJLRGSxiHwqIpEiUltE5orIXyLymYgU8ztOk8eouuTQti0UK+baTDz6qA1AZEyYpNuCW0R6p7P4tOweXESqAb2BBqp62GvTcTNwGfCaqo4RkbeBO4Fh2T2eKSD27YNu3VydxLXXuoGJypTxOypjCrSMuvt4HdA0lkk6yzIbQwkROQ6UBLYA7YFbveUjgf5YsjAAK1fCVVe5u55eeQUeftjqJozJBRkli+nkTEIISlU3icjLwHpcw7/vgQXAHlWN91bbCAStUhWRHkAPgJrWnXTB97//wc03u7Gvp0yBdu38jsiYQiOjZNFJVY+F6+AiUg64CqiN61Lkc+DSIKsGTViqOhwYDhAXFxe2pGZSmjDBh4P+979uaNNzzoGvv4aYGB+CMKbwyqg2cJeIfCUid4lIOG5YvwhYo6rbVfU48CVwHlBWRJISWXVgcxiObbIoV9tZxMe7XmJ793adUs2caYnCGB9klCzuAHYCA4CNIvKbiAwUkZwaY3I9cK6IlBQRATrgOiz8CbjeW6cLMD6HjmdyQFxudfBy6JCrwB4yBB55xFVoR0Xl0sGNMYHSTRaq+qWq3qWq1YAWwFdAJ2CWiGwVkREicr2IRGfl4Ko6FxgH/AYs8uIZDjwBPCIiq4AKwPtZ2b8Jj825cZ23cyd06ODGcB0yxFVmF7Ge8Y3xi6hmvqhfRCrh6hYuwyWP11R1YA7HlilxcXE6f/58P0MoNPr3d1PYbNwIHTvCmjUwerQbyc4YExYiskBVMywvyFKySHWgIkB5Vd2erR1lkyWL3CPi2sSFxapVcNFFrguPSZOgTZswHcgYA6Eni5CGVRWRL9NYdBRYAXyYidiMCW7JEpcojh93XYrnyR4LjSmcQh2D+2rc7auBrZ+SXivwsIi0VdWFORyfKSwWLYL27V0biunToUEDvyMyxgQItSOdIUAiMAEYhLs7SYF3gSlANOBrnYXJx/74Ay68EIoXh59/tkRhTB4U6pVFHeArVb0xaYbXj1NVVb3YK6bKqdtpTWGyeLG766lECVf0dMYZfkdkjAki1CuLC4EYEYkEEJHiQE1vPsDvQPmcD8/kRTk2ntCKFa6OonhxmDbNEoUxeVioVxaLgObAVhHZgGtVHQ3M9ZbH4fpwMoVAjrSzWLPGXVEkJrpEcfrpObBTY0y4hHpl0Q1Yg0sQDYDSwN9AdxFJev58WCI0eU6221hs3eraURw65DoErFcvJ8IyxoRRSFcWqrpMRM4CzsX1ALsJ+EVVE7xVHglTfKag2bcPLr0UtmxxiSI21u+IjDEhCLUYClVNEJHlwAZvVjURQVXXhyc0k1dl+cri6FG4+mp3m+yECW68bGNMvhBqo7yOuP6ZUo8roaHuwxQcVatmod4iMdGNbvfTT/Dxx+7qwhiTb4T6Rf82rlI7NRuirBDasiULGz3zDHz6Kfz73/Cvf+V4TMaY8Aq1grsC8B1QWlVPCZzCGJspKN5/H158Ee6+G554wu9ojDFZEOqX/X9x7SqqeeNOGBOan3+Gnj3h4oth6FAbL9uYfCrUZPEFUBU3MFG8iCR4U3wG2xVqTz/9NBUrVuS0004D4KuvvqJGjRpERUXx+++/+xzdybZu3UqbNm2Ijo7m0UcfzdI+Ro8eTadOndyLNWtc9+JnnAGffQZFc796a8aMGZx11lm5flxjChxVzXACluD6hjppCmX73JiaNWumua1WrVoaGRmppUqVSp7uv/9+VVVdv369RkZG6tatW5PXr1Onjn799dfZPi6gf/31V7b3k9rAgQP1mmuu0cTExKDLu3TpohEREQpRGhUVpWeffbb27dtX9+zZc/LK+/apNmyoWq6c6sqVOR6rMSZnAPM1hO/YUH/q1QTmAX2APTmdsPKziRMnctFFF500f926dVSoUIFKlSqlmHf22WfnZniZsm7dOho0aEB6JY19+vShaNHn6dv3CIsWLaJPnz60bt2auXPnUqpUKbeSqrvzaelS+O47qFs3rHEnJCRQxEbRMyasQi2GGg4cA+ao6h+BUxhjy7emTJlCx44d2bx5M1FRUdxyyy1ERUWRkJBAo0aNON3r2mLz5s1cd911nHrqqdSuXZs333wzeR8JCQm8+OKLnH766URHR9OsWTM2bNhAG28woEaNGhEVFcVnn33Gjh076Ny5M2XLlqV8+fJccMEFJCYmBo1t9uzZNG/enDJlytC8eXNmz54NQNeuXRk5ciSDBw8mKiqKKVOmpHl+/ftDZGQkzZs3Z8KECezcuZMPP3RDmowYMYLz69SBL75AX3qJh7/5hkqVKlGmTBliY2NZvHhx8vF69uxJx44diY6Opm3btqxbty75GMuXL6djx46UL1+es846i7FjxyYv69q1K/feey+XXXYZpUqV4qeffuLbb7+lQYMGREdHU61aNV5++WUApk2bRvXqJ27kW7ZsGe3ataNs2bKcffbZTJgwIcV+77//fi6//HKio6Np2bIlf//9d8Z/cGMKg1AuP4CFQDywF/gDN2b2b8CCULbPjcmvYqgffvgh6LKffvpJq1WrlmIeAcVHCQkJ2rRpUx0wYIAePXpU//77b61du7b+73//U1XVwYMHa8OGDXX58uWamJioCxcu1B07dpy0H1XVvn376j333KPHjh3TY8eO6fTp04MWJe3cuVPLli2rH330kR4/flw/+eQTLVu2bPJ+u3Tpok899VSa55u0vEqVlPNvv/12vfHGG1VV9cPHHtPWoHrjjfq/yZO1adOmunv3bk1MTNSlS5fq5s2bk/cVFRWlP//8sx45ckR79+6trVu3VlXVAwcOaPXq1fWDDz7Q48eP64IFC7RChQq6ePHi5G1Lly6tM2fO1ISEBD18+LCedtppOn36dFVV3bVrly5YsOCkv8OxY8f09NNP1xdeeEGPHj2qU6dO1aioKF2+fHnyfsuVK6dz587V48eP66233qo33XRTmu+HMQUBIRZDhXplEYu7CokGzgEaB0yF2tVXX03ZsmWTp3fffTek7ebNm8f27dt59tlnKVasGHXq1OHuu+9mzJgxALz33ns8//zznHXWWYgIjRo1okKFCkH3FRERwZYtW1i3bh0RERFccMEFQYuSvvnmG+rWrcvtt99O0aJFueWWW6hXrx4TJ07M1DmnHr22atWq7Nq1CzZtgrffhpIl4f33iShWjP3797N8+XJUlfr161MloMvayy+/nDZt2lC8eHFeeOEF5syZw4YNG5g0aRIxMTF069aNokWL0rRpU6677jrGjRuXvO1VV11F69atOeWUU4iMjCQiIoKlS5eyb98+ypUrR9OmTU+K+5dffuHAgQP07duXYsWK0b59ezp37synn36avM61115LixYtKFq0KLfddhsLF9p4XsZA5joSDDZ1D1Nc+cbXX3/Nnj17kqe77747pO3WrVvH5s2bUySaF198ka1btwKwYcOG5OKqjDz++OOcccYZdOrUiTp16jBo0KCg623evJlatWqlmFerVi02bdoU0nGSLFiQ8vWmTZsoX64c3HILHDvmOgaMiqJ9+/b06tWL+++/n8qVK9OjRw/27duXvF2NGjWSn0dFRVG+fHk2b97MunXrmDt3bor3ZvTo0fzzzz9BtwX44osv+Pbbb6lVqxZt27Zlzpw5Qc+/Ro0anHLKiY996vNPunMNoGTJkhw4cCBT740xBVVIyUJVR6Y1hTvAgqpGjRrUrl07RaLZv38/3377bfLyUMvLo6OjeeWVV1i9ejUTJ07k1VdfZerUqSetV7Vq1RT1AgDr16+nWrXUvbik78orTzw/cOAAU6ZM4YK9e2HGDOja1Q1k5OnduzcLFixgyZIlrFy5kv/85z/JyzZs2JBiP7t27aJq1arUqFGDtm3bpnhvDhw4wLBhw5LXT33l1Lx5c8aPH8+2bdu4+uqrufHGG0mtatWqbNiwIUV9TlbO35jCKM1kISLdRCTdu6VEpKiIFPqri6xo0aIFpUuX5qWXXuLw4cMkJCSwePFi5s2bB8Bdd93FM888w19//YWq8ueff7Jz504AKleuzOrVq5P3NWnSJFatWoWqUrp0aYoUKRL07qDLLruMlStX8sknnxAfH89nn33G0qVL6dy5c6bjP3r0KAsWLODqq6+mXLFidPv+e7jrrhSdA86bN4+5c+dy/PhxSpUqRWRkZIq4vv32W2bOnMmxY8d45plnaNmyJTVq1KBz586sXLmSjz/+mOPHj3P8+HHmzZvHsmXLgsZy7NgxRo8ezd69e4mIiEh+D1Jr2bIlpUqVYvDgwRw/fpxp06YxceJEbr755kyfvzGFTXpXFu8Dm0RkqIjcLCKNRKSO93iziAzDdVUeWiF9AXXFFVcQFRWVPF1zzTUhbVekSBEmTpzIwoULqV27NhUrVuSuu+5i7969ADzyyCPceOONdOrUidKlS3PnnXdy+PBhAPr370+XLl0oW7YsY8eO5a+//uKiiy4iKiqKVq1acd9999GuXbuTjlmhQgUmTZrEK6+8QoUKFRg8eDCTJk2iYsWKIZ/v4MGDgWjKly/PHXfcQbP69Zl95AilGjSAN95Ise6+ffu4++67KVeuHLVq1aJChQo89thjyctvvfVWBgwYQPny5VmwYAGjR48G3JXS999/z5gxY6hatSqnnXYaTzzxBEePHk0zro8//piYmBhKly7N22+/zahRo05ap1ixYkyYMIHJkydTsWJF7rvvPj766CPq2XgaxmRIXGV4kAXuiqE/rgPBYCsJbnS8fqr6YbgCDFVcXJzOT13zasJCxDWlQNWVSf3wA/z6a6bGpujatSvVq1fn+edtzCxj/CQiC1Q1LqP10ixmUtUPRGQkcAVwOe6OqHK4Rnl/ApOAiXpiACRT2AwbBpMmweuv2yBGxhRw6dZJeInga28y5oQVK+DRR+GSS6B3b7+jMcaEmQ1cZDJt07p4uLGLu+vpgw+y1JPsiBEjcj4wY0zY2HgUOezIkSMcPHjQ7zDCasFTX8Lcua7L8YBGdgXFrl270uwuxZjCypJFDjl8+DCvvPoaVarX5NE+ff0OJ3z+/JPho0vBjTdCAb3l9IILO3Bmg1i++OILSxrGeHwvhhKRssB7QEMWXieZAAAgAElEQVTcXVfdgRXAZ0AMsBa4UVV3+xRiug4ePEi//gN49/0PKXLamSTUasmWLVtYtGiR36HlvPh4uO02Xiy3lUW9voaCeI7Anr37OVS7Dd0fepKHH+9Ln0cepFevXn6HZYyv0rx1NtcCcHdczVDV90SkGFASeBLYpaqDRKQvUE5V0x2P069bZ9955x163nsfJWo3pcLlD3NoxSz0t3FUO61SxhvnN9u3w9atrIusR60zfP+dETYr12ygwk0vIFKEbV8MQA/uYs3fq07qYsSYgiDUW2dD7XU2FvgZ2AckBEzxoWyfzn5LA2vwklbA/BVAFe95FWBFRvvyo9fZJOvWrdNud/XQktFlNLLqWXpH97t8iyVsli1TLV5c9brrFPwOJrxqnVFPS9dupNHlKujA55/Xffv2+R2SMWFDiL3OhnRlISJ/4HqbDZZsslzvISKNcWNlLAUaAQuAB4FNqlo2YL3dqlouyPY9gB4ANWvWbJa636Pctn79evo/9wLNmjbh/nt7+hpLjkpMhDZt3GBGS5ciVU7D5wvSsLq7533UrFGNh3r3Jjo62u9wjAmrUK8sQk0WB4HVwAO4RnnJG2k2BkASkTjgF6C1qs4VkTdwVy8PhJIsAlkL7jAaPhzuuQc+/BC6dj3RgtsYk+9luwV3KlOB46o6LVtRnWwjsFFV53qvxwF9ga0iUkVVt4hIFWBbDh/XhGrrVnjiCWjXDrp08TsaY4xPQk0Wm4AeIjKBE6PmAaCqA7N6cFX9R0Q2iMhZqroC6IArkloKdAEGeY/js3oMk02PPAKHDrmuPbLQ+M4YUzCEmizu8R474/qJAteRoAJZThaeB4DR3p1Qq3GDKp0CjBWRO4H1wA3ZPIbJih9+gE8+gWefdQMaebLQo7kxJp8Ltc5iBMF7nkVVu+VwTFlidRY57Ngx1zlgQoJrTxEZ6XdExpgwyNE6C1Xtmu2ITP7y+uuus8Bvvz0pUVxxBWRy2G5jTD4X8m2vItJKREaJyAwR+VhEWmW8lcmXNm6EgQPhqqvg0ktPWtyjhw8xGWN8FdKVhYhcjBu/ImmsytbALSJyhapODldwxiePP+6Kn157LejiZs1yOR5jjO9CvbLoj6uzeB3o6T0mAs+GJyzjmxkzYMwYd7ts7dpBV6lWLZdjMsb4LtS7oRoAn6rqI0kzRKQCcHVYojL+SEyEhx6C6tWhTx+/ozHG5CGhJos9wNkiEqmqR0QkEjjbm28KipEj4bffYPRoKFnS72iMMXlIqMVQPwBNcS2rFwNbgSbefFMQ7N8PTz4J554Lt9zidzTGmDwm1CuLvrgriZa4IimAud58UxC89BL88w989ZW11DbGnCTUdhY7ROQ8II4TAxLN11Ba9Jm8b9MmePVVN/Lduef6HY0xJg9KM1mISE1gn6ru8Z6DK37a6j2vISKo6vpwB2nCrF8/Nwreiy+GtPqECWGOxxiT56R3ZbEGd4vso7griWBXEZrBPkxet2SJ63r8wQfTvFU2NWtnYUzhk14Ft3hT6teBU5YHPjJ5RN++EB0NTz0V8iZxGQ/AaIwpYNK8KggcAS87o+GZPGzGDJg0CQYNggoVQt5s8+YwxmSMyZNCSgIi8qOI3Jtq3jXeyHYmP1J1t8pWrQq9e2dq0/79wxOSMSbvCvWKoR1QN9W8C4FeORqNyT2TJ8PMmfDMM1CiRKY2HTAgTDEZY/KsdCunReSDgJedAl6fAlwGHA5XYCaMEhNdHUWdOtC9u9/RGGPygYzuZOqKu+NJcY3xGqRa/r8wxGTC7fPPYeFCGDUKihXzOxpjTD6QUbJIKnDoh2uxnZQcEoANwOdhisuES0KCq3Q4+2zXCM8YY0KQbrJQ1QEA4rp/+EVVv8uNoEwYjRkDy5fDuHFQpEjG6xtjDKF39zFARGqKyM3AqQS0v1DVN8MVnMlh8fGudjo2Fq65Jsu7qVIlB2MyxuQLoY6Udx3wMVA8yGJLFvnFJ5/AX3/Bl1/CKVlvOmPtLIwpfEL9xngGiMD1CyXAMiAe+DlMcZmcFh8Pzz0HTZrA1dkbs8raWRhT+ITar1M94AtgM/CgqjYUkR+A+WGLzOSsTz6BVavg66+tC3JjTKaFemURD2wHDgCISG3gEG48bpPXJSS4HmUbNYIrr8z27uzKwpjCJ9Qri81AFeBXXDHUCqAIJ7orN3nZuHGwYgWMHZsjVxVVq1q9hTGFTahXFsOBLbiiqCW4JKPAwDDFZXJKYiI8/zzUrw/XXZcju9yyJUd2Y4zJRzK8shDXyGIscFRVD4lIU1xL7h2quincAZpsmjABFi+Gjz/O1h1QxpjCLdRvj1XASwCqelxV/7BEkQ+ourqK00+31trGmGzJ8MpCVVVEfgfK5UI8Jif9+CPMmwfvvANFbUBDY0zWhfoN8hPwmIh8CMwGjiYtUNWPwhGYyQGDBrnm1l26+B2JMSafCzVZ9MFVaN/hTYEsWeRF8+fDlCkweDAUD9bwPuv69cvR3Rlj8oFQk8V0XLIICxEpgmvgt0lVO3vtOMYA5YHfgNtV9Vi4jl8gDRoEZcvCPffk+K6tnYUxhU+oHQm2C3McD+K6ECntvX4JeE1Vx4jI28CdwLAwx1BwrFjh+n968kkoXTrj9TPJ2lkYU/iEOgZ3goi8kmpeHxH5I7sBiEh14HLgPe+1AO2Bcd4qI4HsdWZU2LzyihvU6IEHwrL7+dbJizGFTqi3zgoB3ZJ76gINcyCG13F1Ione6wrAHlWN915vBKoFDUqkh4jMF5H527dvz4FQCoCtW+Gjj1ylduXKYTnEggVh2a0xJg9LN1mIyGoRWe297J70WkTWAt2APdk5uIh0BrapauDXT7D+KILWl6jqcFWNU9W4U089NTuhFBxDhsCxY/Doo2E7RA50L2WMyWcyqrOI8R4VV5+QugD8vWwevzVwpYhcBkR6+38dKCsiRb2ri+q4vqlMRg4ehKFD4aqr4Mwz/Y7GGFOAZFQMdSGu/kBwdQgXelMboLaqPpGdg6vq/6lqdVWNAW4GflTV23DtOq73VusCjM/OcQqNDz+EXbvg8cf9jsQYU8BkNAb3zwAiciGwUVX/zpWo4AlgjIg8D/wOvJ9Lx82/EhLgtdegVSs47zy/ozHGFDChtrOYBdwuIg8CUQHzVVXvzIlAVHUaMM17vhpokRP7LTQmTIDVq+Gll/yOxBhTAIWaLD4EbvWeB1ZAK64NhPHba69BTEy2h0w1xphgQk0WV+H6gxqLuwMqbK25TRbMnw8zZsCrr+ZKh4GbrL9hYwqdUL9ZtgCzVLV7OIMxWfTaaxAdDXfmzkXeggWuFbcxpvAItVHem8AVInKtiNQRkZpJUziDMyHYtMkNl3rnnWHp2iOY4cNz5TDGmDwk1CuL/+KKnj5PNV8zsQ8TDm+95YZO7d071w45cWKuHcoYk0dkZpxNCTLZOJ1+OnzY/cy/8kqoXTvXDnvFFbl2KGNMHhFqr7OWFPKiTz6BnTvhwQdz9bCTJuXq4YwxeUDISUBEIkSkk4jcKyLFvTqLkuEMzqRDFd54A2JjoW1bv6MxxhRwoXZRXgNYCEzG1V+UBVYBz4UvNJOuadNg0SJ3VSHB+l40xpicE+qVxWtAfWAHIKq6FTd63iXhCsxk4M03oWJFuPXWjNc1xphsCjVZtAEmAZ8EzFsF2K2zfli71nXv0aMHREb6HY0xphAINVkokHoM7BjgQI5GY0IzbJgreurZ05fDd+7sy2GNMT4KtY3EfKAz7moCEfkc6Ii72jC56fBheO891wdUjRq+hGDtLIwpfEK9sugDHAIa4NpXXAfsBZ4KU1wmLZ984sasCNP42qGwdhbGFD6htrNYIiL1gdtxxU9rgdGquiV8oZmTqLphU885B9q08S2MHj18O7QxxichJQsRiQQOA6+oqnrzSotIpKoeCWeAJsCsWbBwIbzzjq+3yzZr5tuhjTE+CbUY6idcvQUAIiLAL8DUcARl0jBkCJQpA7fd5msY1ar5enhjjA9CTRYNgZ+Triq8xxlAbLgCM6ls2QJffAHdu0OpUn5HY4wpZEJNFoeBM1LNO9Obb3LDu+9CfDzce6/fkRhjCqFQb51dCHQQkS+BuUBLXEO9KeEKzAQ4ftzVU1x8MdSt63c0xphCKNRk8SxwAXA1bohVwQ2z+kyY4jKBxo+HzZvh7bf9jsQYU0iFeuvsLyISB9yLu3V2DfC2qi4JY2wmyVtvQUwMXHaZ35EYYwqpkEe58xJDrzDGYoJZssT1MDtoEBQp4nc0gOuWyhhTuITazqIM8AjQBIgKWKSq2iEcgRnP0KFQvLgbYzuPsHYWxhQ+oV5ZfILrjjx1SzDN2XBMCvv3w0cfwU03ue7I84i4OFeFYowpPEJNFm2BXcBbwB4sSeSOjz+GAwfg/vv9jiQFSxTGFD6hJouVwB+q2j+MsZhAqq5iu1kzaN7c72hS6N/fTcaYwiPURnmfATeIyCMi0l5E2iRN4QyuUJs+HZYudVcVeWzY1AED/I7AGJPbQr2y+Deu6Ok/qeZrJvZhMuOtt6B8ebj5Zr8jMcaYkL/o12P1FLln0yb48kt4+GEoUcLvaIwxJuRGeTHhOLiI1AA+Ak4DEoHhqvqGiJTHFX3F4MbOuFFVd4cjhjxp+HBITLR+oIwxeUaodRaISBkRuUtEnhORO0WkbA4cPx54VFXrA+cC94tIA6AvMFVV6+K6Qe+bA8fKH44dc8ni0kuhTh2/ozHGGCD0Rnln4sa0OC1g9kARaa+qK7J6cG+kvS3e8/0isgyohut/qp232khgGvBEVo+Tr3z1FfzzT567XTZQlSp+R2CMyW2hXlm8DFQB/sAVDy30Xg/OqUBEJAbXQnwuUDlpyFbvsVIa2/QQkfkiMn/79u05FYq/hgxxVxSXXOJ3JGmydhbGFD6hJotzge9Vtamq3qqqzYDvgFY5EYSIRAFfAA+p6r5Qt1PV4aoap6pxp556ak6E4q8//oCZM+G+++CUkEsIc521sTCm8An1G0mAY6nmHePk7j8yTUQicIlitKp+6c3eKiJVvOVVgG3ZPU6+MGSIu/upe3e/IzHGmBRCTRYLgMtFZKqIDBGRKUBnAsblzgpvLO/3gWWq+mrAoglAF+95F2B8do6TL+zaBaNHw7/+BeXK+R1NuuzKwpjCJ9R2Fn2An4ELcRXPAuwl+3cptQZuBxaJyEJv3pPAIGCsiNyJa+NxQzaPk/e9/z4cPgwPPOB3JBmqWtXqLYwpbEJtZ/GniNTHfbHH4No+jEqqhM4qVZ1J2kVZhafr84QE1xV527Zwzjl+R5OhLdn6qxtj8qN0k4WIFAVKAodU9R+87j6S5otIUVWND3+YBdw338DatfDyy35HYowxQWVUZ/EqsB2okWp+VWAr8Eo4gip03ngDatSAq67yOxJjjAkqo2RxGTBNVdcEzlTV9biW1TYodHb9+Sf8+KOrqyhqfTIaY/KmjJJFNVz9RDCbgOo5Gk1h9MYbULIk3HWX35EYY0yaMkoWB4GTRt7xbnlt4S03WbVtm7tdtmvXPH+7bKB+/fyOwBiT2zJKFnOARiIyQkTqikgxEakLjABiveUmq955B44ehd69/Y4kU6ydhTGFT0aF5IOAS3C3zN6ealki8FI4gioUjh51t8tedhmcdZbf0WSKtbPIu44fP87GjRs5cuSI36GYPCYyMpLq1asTERGRpe3TTRaqOktEbgeGAOUDFu0GHvDaSZisGDXK9S778MN+R5Jp87PVbt+E08aNG4mOjiYmJgbJY8PxGv+oKjt37mTjxo3Url07S/vIsLsPVR2Du3W2E/Av77GGqn6SpSMaN7DRf/4DTZpAh/zX9nDBAr8jMGk5cuQIFSpUsERhUhARKlSokK0rzlBbcB8GpmT5KCaliRNhxQr49FPIh//UV14JaoPs5lmWKEww2f1c5N1+sAuywYMhJgauv97vSIwxJiSWLHLbrFkwezY88og1wjMFUpEiRWjcuDGNGjWiadOmzJ49O8f2HRMTw44dO1LM27NnD0OHDk1+PW3aNDp37pxjxwxVVFRUrh8zN1myyG0vvgjly9uYFabAKlGiBAsXLuSPP/7g3//+N//3f/930joJCQk5drzUySJUORlDboiP97cbPvtpm5vmzYNvv4UXXoBSpfyOxhR0Dz0ECxdmvF5mNG4Mr78e8ur79u2jnNfgdNq0aQwYMIAqVaqwcOFCli5dyqhRo3jzzTc5duwYLVu2ZOjQoRQpUoR7772XefPmcfjwYa6//noGDBiQYr+HDx/mmmuu4brrrmPq1Kn8/fffNG7cmI4dO3L55Zdz4MABrr/+ehYvXkyzZs0YNWoUIkJMTAzdu3fn+++/p1evXtSrV4+ePXty6NAhTj/9dD744APKlStHu3btePnll4mLi2PHjh3ExcWxdu1aDh06RNeuXVm+fDn169dn7dq1vPXWW8TFxQHw1FNPMWnSJEqUKMH48eOpXLlyirgPHDjAAw88wPz58xER+vXrx3XXXUdUVBQHDhwAYNy4cUyaNIkRI0bQtWtXypcvz++//07jxo356quvWLhwIWXLlgXgjDPOYNasWZxyyin07NmT9evXA/D666/TunXrrP2N02DJIjcNHOhaavfq5XckxoTN4cOHady4MUeOHGHLli38+OOPyct+/fVXFi9eTO3atVm2bBmfffYZs2bNIiIigvvuu4/Ro0dzxx138MILL1C+fHkSEhLo0KEDf/75J7GxsYD7wr355pu54447uOOOO+jYsSOLFy9moZcYp02bxu+//86SJUuoWrUqrVu3ZtasWZx//vmAa28wc6a76z82Npb//ve/tG3blmeffZYBAwbwejrJcOjQoZQrV44///yTxYsX07hx4+RlBw8e5Nxzz+WFF16gT58+vPvuuzz99NMptn/uuecoU6YMixYtAmD37t0Zvp8rV65kypQpFClShMTERL766iu6devG3LlziYmJoXLlytx66608/PDDnH/++axfv56LL76YZcuWhfLnCpkli9zy228waRI89xyULu13NNmyaZPfEZiQZOIKICclFUMBzJkzhzvuuIPFixcD0KJFi+T7/KdOncqCBQto3tz1KHT48GEqVaoEwNixYxk+fDjx8fFs2bKFpUuXJieLq666ij59+nDbbbelGUOLFi2oXt11Xde4cWPWrl2bnCxuuukmAPbu3cuePXto27YtAF26dOGGG9IfZ23mzJk8+OCDADRs2DA5JoBixYol15U0a9aMH3744aTtp0yZwpgxY5Jflwuhm58bbriBIkWKJMc+cOBAunXrxpgxY5LPZcqUKSxdujR5m3379rF//36io6Mz3H+oLFnkloEDoWzZfDESXkYWLHCtuI3JSKtWrdixYwfbt28HoFRA8auq0qVLF/7973+n2GbNmjW8/PLLzJs3j3LlytG1a9cU7QNat27N5MmTufXWW9O8HbR48eLJz4sUKZKivL9UCEXARYsWJTExESDFsTWde8YjIiKS40l9zMDtg8UcOC91W4jAeFu1asWqVavYvn07X3/9dfKVS2JiInPmzKFEiRIZnltWWQV3bpg3D8aPd621y5TxO5psGz7c7whMfrF8+XISEhKoUKHCScs6dOjAuHHj2LZtGwC7du1i3bp17Nu3j1KlSlGmTBm2bt3K5MmTU2w3cOBAKlSowH333QdAdHQ0+/fvz3RsZcqUoVy5csyYMQOAjz/+OPkqIyYmhgVe69Nx48Ylb3P++eczduxYAJYuXZpcnBSqTp06MWTIkOTXScVQlStXZtmyZcnFTGkREa655hoeeeQR6tevn/y+pt7vwpyuq8KSRfipwmOPQaVK+bJrj2AmTvQ7ApOXJdVZNG7cmJtuuomRI0cmF6MEatCgAc8//zydOnUiNjaWjh07smXLFho1akSTJk04++yz6d69e9CK2tdff50jR47Qp08fKlSoQOvWrWnYsCGPP/54pmIdOXIkjz/+OLGxsSxcuJBnn30WgMcee4xhw4Zx3nnnpbhV97777mP79u3Exsby0ksvERsbS5lM/AB8+umn2b17Nw0bNqRRo0b89NNPAAwaNIjOnTvTvn17qlSpku4+brrpJkaNGpVcBAXw5ptvMn/+fGJjY2nQoAFvv/12Zt6GkEh6l1X5SVxcnM7Pi50WjR8PV1/tOg28916/o8kRV1xhCSOvWrZsGfXr1/c7jAIrISGB48ePExkZyd9//02HDh1YuXIlxYoV8zu0kAT7fIjIAlWNy2hbq7MIp+PH4YknoF69AjW40aRJfkdgjD8OHTrEhRdeyPHjx1FVhg0blm8SRXZZsgind991fUCNHw9Z7BbYGJN3REdHkydLMHKB1VmEy+bN8OSTcOGFrtzGGGPyMUsW4dKrlxvg6J138mXPssYYE8iKocLhiy/gq69g0CCoW9fvaIwxJtvsyiKn7djhriqaNIFHH/U7mrDwoUNPY4zPLFnkpPh4uOkm2L0bPvigwHZBbrfNmvQkdVHesGFDbrjhBg4dOgTkbBfet9xyC7Gxsbz22mt07do1RcM5Ex6WLHJSnz7w44+uniKgg7GCxurrTXqS+oZavHgxxYoVy/EGYv/88w+zZ8/mzz//5OEC0tA1P7BkkVNGjoTXXnN9P3Xp4nc0YdWjh98RmPziggsuYNWqVSnmpR6cqFevXowYMQKAb7/9lnr16nH++efTu3fvoIMYderUiW3bttG4cePkrjqSTJ06lSZNmnDOOefQvXt3jh49yq+//sq1114LwPjx4ylRogTHjh3jyJEj1KlTJ4fPuOCyZJEThg6Fbt2gfXt45RW/owm7Zs38jsCEqn9/dzNe6impI8j+/d0Ebl6wdQOXb94c+rHj4+OZPHky55xzTkjrHzlyhHvuuYfJkyczc+bM5M4HU5swYQKnn346Cxcu5IILLkixfdeuXfnss89YtGgR8fHxDBs2jKZNm/L7778DMGPGDBo2bMi8efOYO3cuLVu2DP2ECjlLFtmRmAj9+sH997ta30mTCkXju2rV/I7AhKp/f9c9Weop6Us/MFls3hx83cDlofQ2nNQ3VFxcHDVr1uTOO+8MKdbly5dTp06d5C7Mb7nllsycKitWrKB27dqceeaZgOtyfPr06RQtWpQzzjiDZcuW8euvv/LII48wffp0ZsyYkSLZmPTl2RpYEbkEeAMoArynqoN8Diml2bNdx4C//uqGSH3nnQJboW1MZgSOZxFMYPffcKJL7uz2U5fe9hdccAGTJ08mIiKCiy66iK5du5KQkMDLL7+crWMWJnnyykJEigBvAZcCDYBbRKSBr0EdOQJLl8Kbb8LFF0Pr1rBxI4wYAe+9Z4nCmBDVqlWLpUuXcvToUfbu3cvUqVMBqFevHqtXr2bt2rUAfPbZZ5nab7169Vi7dm1yHUlgl+Nt2rTh9ddfp1WrVpx66qns3LmT5cuXc/bZZ+fciRVwefUbrgWwSlVXA4jIGOAqYGm6W2XFBx/Af/4DCQluSqLqipkSEmDvXgjsL79uXVf89NhjkIO3AxpTGNSoUYMbb7yR2NhY6tatS5MmTQB3RTJ06FAuueQSKlasSIsWLTK138jISD788ENuuOEG4uPjad68OT179gSgZcuWbN26lTZt2gBuONVKlSqlOXiSCUJV89wEXI8rekp6fTswJMh6PYD5wPyaNWtqlowfr/0ajA1SUqtapcQu1e7dtV/Lydqv/XTVDz/UKpWOB123Xz+3uypVVDdtUp0wIVjpr5uSlnfu7Lbp3Dn4eoHLJ0xw26W1z6TlVaq4bfr1C75e4PLAmDNzTuXKZe2tNuG3dOlSv0PIlv3796uqamJiot5777366quv+hxRwRLs8wHM1xC+l/PkeBYicgNwsare5b2+HWihqmmOSZpnx7MwJhfl9/EsXnvtNUaOHMmxY8do0qQJ7777LiVLlvQ7rAKjII5nsRGoEfC6OpCJm/aMMfnRww8/bA3t8qg8WcENzAPqikhtESkG3AxM8DkmY/KFvFhaYPyX3c9FnkwWqhoP9AK+A5YBY1V1ib9RGZP3RUZGsnPnTksYJgVVZefOnURGRmZ5H3m1GApV/Rb41u84jMlPqlevzsaNG9Ns/WwKr8jISKpXr57l7fNssjDGZF5ERERyC2hjclKeLIYyxhiTt1iyMMYYkyFLFsYYYzKUJxvlZYWIbAfWZXHzisCOHAzHT3YueU9BOQ+wc8mrsnMutVT11IxWKjDJIjtEZH4oLRjzAzuXvKegnAfYueRVuXEuVgxljDEmQ5YsjDHGZMiShTPc7wBykJ1L3lNQzgPsXPKqsJ+L1VkYY4zJkF1ZGGOMyZAlC2OMMRkq9MlCRC4RkRUiskpE+vodT1aJSA0R+UlElonIEhF50O+YskNEiojI7yIyye9YskNEyorIOBFZ7v1tWvkdU1aJyMPeZ2uxiHwqIlnvwjSXicgHIrJNRBYHzCsvIj+IyF/eYzk/YwxFGufxH+/z9aeIfCUiZcNx7EKdLESkCPAWcCnQALhFRBr4G1WWxQOPqmp94Fzg/nx8LgAP4rqnz+/eAP6nqvWARuTTcxKRakBvIE5VGwJFcOPM5BcjgEtSzesLTFXVusBU73VeN4KTz+MHoKGqxgIrgf8Lx4ELdbIAWgCrVHW1qh4DxgBX+RxTlqjqFlX9zXu+H/elVM3fqLJGRKoDlwPv+R1LdohIaaAN8D6Aqh5T1T3+RpUtRYESIlIUKEk+Gr1SVacDu1LNvgoY6T0fCVydq0FlQbDzUNXvvTGAAH7BjSya4wp7sqgGbAh4vZF8+gUbSERigCbAXH8jybLXgT5Aot+BZFMdYDvwoVek9p6IlPI7qKxQ1U3Ay8B6YAuwV1W/9zeqbKusqlvA/dgCKvkcT07oDkwOx44Le7KQIPPy9b3EIhIFfAE8pKr7/I4ns0SkM7BNVRf4HUsOKAo0BYapahPgIPmjqOMkXiuV9IEAAARiSURBVHn+VUBtoCpQSkT+5W9UJpCIPIUrjh4djv0X9mSxEagR8Lo6+ejSOjURicAlitGq+qXf8WRRa+BKEVmLKxZsLyKj/A0pyzYCG1U16QpvHC555EcXAWtUdbuqHge+BM7zOabs2ioiVQC8x20+x5NlItIF6AzcpmFqPFfYk8U8oK6I1BaRYrgKuwk+x5QlIiK4svFlqvqq3/Fklar+n6pWV9UY3N/jR1XNl79gVfUfYIOInOXN6gAs9TGk7FgPnCsiJb3PWgfyaWV9gAlAF+95F2C8j7FkmYhcAjwBXKmqh8J1nEKdLLxKoV7Ad7gP/lhVXeJvVFnWGrgd90t8oTdd5ndQhgeA0SLyJ9AYeNHneLLEuzoaB/wGLMJ9d+Sb7jJE5FNgDnCWiGwUkTuBQUBHEfkL6Oi9ztPSOI8hQDTwg/d//3ZYjm3dfRhjjMlIob6yMMYYExpLFsYYYzJkycIYY0yGLFkYY4zJkCULY4wxGbJkYYxHREaIiIpI/1w+7oXecZ/zXrfzXq/Nwr7u9rbtmtNxmsLNbp01hYL3xVsrnVUuxHVj0QLXS+z/ciMuABGZA8QB1VV1q4i0A34C1nmNEzOzr+K4luOHgDqqmpDD4ZpCqqjfARiTSz4AynvP7wWK4bpG2ejN26iq04BPcjMoEWmM61L+e1Xdms19FVXVoyLyFXA3ruv9fD0eiMk7rBjKFAqqOlBVH1LVh4DD3uwhSfNUdVXqYigR6e+9niwio0TkkIj8KiJniMhwETngDQSU3N+TiNQUkTEisklE9ojI9yLSMJ3QOnuPPwdbKCKPiMg/3oA3jwfMT4r1HW/gnmPA+an21fmkHRqTRZYsjMnYxUAFXP9IzYH5uC7g/wTOBt4EEJGSwI/Ajd6yH4B2wE8iUjGNfcd6j8H6WaoJ3AXMBk4FXpL/b++OVasIogAM/0dIp76AEIQoYhXEB0gIVnYB+4CIT5E27xCwSCPkCZLOJhCxkYCVQSJGEGIpYiEET4rZgcuyy0TuQiL+XzOXnZlz9zZ7LnPnzom43xvzElgAXgP1lOEa69FlP6DU4jKU1HYCPKUcNrdDKfzzhPIwr2c+QSnYtAR8A467a1+7a8+AoTN7ainPnwN9f4C1zDyLiNPu/ZaBTzNjDjJztTevJo1rXyZU/w6ThdT2MTMzImqVu++Z+SMi6gO+FjS627V3KGVhZ90biV1j3hroO+tOrq3jFoGbvTFvB+bd7sWW5uYylNTW31E0tsPoS9e+B25kZmRmUL7hb43M+dC1Dwf6zmdej21b/D1wrcY6Gpkj/TWThTSdfeAz8Bg4jIjtiNinFNRaHplTdyutTHgfNdbehDH1nzNZSBPJzF+UokC7lCWjDeAB5cfn45E5R8A7Sh2SuWtAd0W81im15U0Wmox/ypOuWESsAW+ArczcnDPWC+AV8Dwzd6a4PwlMFpKkS3AZSpLUZLKQJDWZLCRJTSYLSVKTyUKS1GSykCQ1mSwkSU0XpGVDNNYP02gAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "%matplotlib inline\n", "\n", "# Define an array time points to estimate the function\n", "times = np.arange(0,12,0.1)\n", "\n", "# Define a function to use to compute the value of the \"ith\" term\n", "# in the series of eigenfunctions that are summed in the solution\n", "def eigenfunction(Pe, B, x, t):\n", " return (B * (B * np.cos(B * x) + Pe/2 * np.sin(B * x)) /\n", " (B**2 + Pe**2/4 + Pe) / (B**2 + Pe**2/4) * np.exp(-1 * B**2 * t)\n", " )\n", "\n", "# Store the results in a list\n", "Cs = []\n", "\n", "# Estimate the concentration for each dimensionless time at x = 1\n", "for t in times:\n", " tau = D * t / L**2\n", " x = 1\n", " \n", " # Get the eigenfunction values for all the eigenvalues\n", " series = eigenfunction(Pe, np.array(betas), x, tau)\n", " \n", " # Sum the series and convert the result to concentration at the point\n", " C = C0 * (1 - 2 * Pe * np.exp(Pe/2 * x - Pe**2/4 * tau) * series.sum())\n", " Cs.append(C)\n", " \n", "# Plot the results\n", "fig, ax = plt.subplots()\n", "ax.set_xlabel('Time (hr)', size = 12, weight = 'bold')\n", "ax.set_ylabel('Concentration (mg/L)', size = 12, weight = 'bold')\n", "ax.set_title('Column Breakthrough Curve', size = 14, weight = 'bold')\n", "ax.plot(times, Cs, ls = '-', c = 'r', label = 'Breakthrough curve')\n", "\n", "# Add a couple of other lines for explanation of behavior\n", "xs = [0, L/U, L/U, 12]\n", "ys = [0, 0, C0, C0]\n", "ax.plot(xs, ys, ls = '-.', lw = 1, c = 'b', label = 'Plug flow')\n", "ax.text(0.5,65,'Effects of Dispersion', fontsize = 12)\n", "arrowprops = {'arrowstyle':'<|-|>'}\n", "ax.annotate('', xy = (5.5,61), xytext = (0.5,61), ha = 'right', va = 'center', arrowprops = arrowprops)\n", "leg = ax.legend()\n", "plt.savefig('breakthrough')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The dispersion in the column caused by differences in flow path lengths acts to broaden the breakthrough curve of the tracer compared with ideal plug flow. \n", "\n", "### Dispersion model fitting\n", "The model is useful for fitting values of dispersion coefficients to observations of the concentration in the effluent. In this case, the model must be run repeatedly at different values of $D$ to meet some fitness criteron, such as minimizing the total squared error." ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.7.0" } }, "nbformat": 4, "nbformat_minor": 2 }